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ABSTRACT 


The  nonlinear  buckling  analysis  of  geometrically  imperfect,  thin  circu- 
lar, cylindrical,  stiffened  shells  under  pure  torsion  and  torsion  combined 
with  axial  compression,  for  various  transverse  and  in-plane  boundary  condi- 
tions, is  first  performed.  A methodology  is  presented  for  predicting  criti- 
cal conditions  (limit  point  loads)  for  such  configurations.  This  methodology 
is  based  on  the  smeared  technique  (closely  spaced  stiffeners),  the  von 
Karman-Donnell  nonlinear  kinematic  relations  in  the  presence  of  initial  im- 
perfections, and  linearly  elastic  behavior.  Hie  computational  procedure 
employs  a Fourier  series  type  of  separated  solution  and  by  employing  the 
Galerkin  procedure  in  the  circumferential  direction  and  identities  of  trigo- 
nometric functions  the  field  equations  are  reduced  to  a system  of  ordinary 
differential  equations.  These  equations  are  then  solved  by  the  finite  dif- 
ference scheme.  Numerical  results  for  numerous  stiffened  and  unstiffened 
configurations  are  presented.  Some  of  these  configurations  are  used  as 
benchmarks  for  the  developed  methodology,  since  numerical  solutions  for  them 
have  been  reported  in  the  open  literature. 

Then,  the  effect  of  initial  geometric  imperfections  on  optimal  stiffen- 
ed configurations  is  assessed.  This  is  accomplished  by  computing  the  criti- 
cal load  at  the  optimum  design  point  as  well  as  at  design  points  in  the  sur- 
rounding space  for  a given  geometric  imperfection.  The  optimum  point  was 
obtained  through  the  use  of  linear  buckling  analysis.  The  comparison  shows 
(in  this  case)  that  the  optimum  point  for  tors  ion- loaded  perfect  stiffened 
cylinders  is  also  the  optimum  design  point  for  same  geometry  imperfect 
stiffened  cylinders  provided  that  the  proper  factor  is  used  with  the  applied 
torsion  to  account  for  imperfection  sensitivity.  Finally,  a sumnary  of  the 

work  performed  under  this  grant  is  presented  at  the  end  of  this  report. 
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1.  INTRODUCTION 

Most  of  the  investigations  reported  in  the  open  literature  (see  Ref.  1 

and  the  cited  references  therein  for  a complete  and  comprehensive  historical 

sketch)  which  consider  buckling  of  imperfect  shell  configurations,  employ  an 

isotropic  constant  thickness  geometry  and  a uniform  axial  compression.  Very 

2 

few  consider  other  constructions  and  even  fewer  torsional  loading.  Loo  , in 

1954,  reported  the  results  of  his  investigation  of  the  effect  of  initial 

imperfections  on  the  critical  condition  for  simply  supported,  constant 

thickness,  isotropic,  thin,  cylindrical  shell  loaded  in  torsion.  His  formu- 

3 

lation  contains  numerous  approximations  and  simplifying  assumptions.  Nash 
removed  some  of  the  simplifying  assumptions  and  extended  Loo's  work  to  the 
case  of  clamped  boundary  conditions.  Both  of  these  investigations  have 
made  important  contributions  to  the  state  of  the  art, but  their  reported 
results  can  only  be  considered  as  qualitative. 

A different  approach,  based  on  Koiter's  initial  postbuckling  theory 
(Ref.  4 and  5), was  employed  by  Budiansky^  in  dealing  with  the  same  problem 
as  Loo  and  Nash.  Results  were  presented  for  classical  simply  supported  and 
two  sets  of  clamped  boundary  conditions. 

The  present  report  is  an  extension  of  the  work  reported  by  the  au- 
thors in  Refs.  1,  7.  These  references  present  the  buckling  analysis  of  im- 
perfect, stiffened,  thin  cylinders  of  finite  length  under  uniform  axial 
compression  and/or  lateral  pressure,  for  axisymmetric  or  at  most  symmetric 
imperfections  and  various  boundary  conditions.  The  present  work  removes 
the  limitation  on  the  imperfection  shape  and  considers  individual  or  com- 
bined application  of  uniform  axial  compression,  lateral  pressure  (which  can 
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be  position  dependent) , and  uniform  end  torsion.  The  methodology  of  Ref.  1 
is  modified  in  order  to  accomodate  the  general  shape  of  the  imperfection  and 
the  new  load  conditions.  Results  of  this  investigation  are  presented  both 
in  graphical  and  tabular  form  for  several  examples. 

In  addition,  the  optimization  of  imperfect  stiffened  cylinders  under 
torsion  is  investigated.  This  is  accomplished  by  first  finding  the  optimum 
geometry  of  the  corresponding  perfect  configuration  by  using  linear  buckling 
analysis  and  a safety  factor  to  account  for  imperfection  sensitivity,  and 
then  by  assessing  the  effect  of  geometric  imperfections  on  this  optimum 
geometry  and  the  surrounding  design  configurations. 
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2.  MATHEMATICAL  FORMULATION 


Given  an  imperfect  (w  ) , stiffened,  thin,  circular  cylindrical  shell 
of  finite  length,  L,  and  various  boundary  conditions  under  the  application 
of  axial  compression,  lateral  pressure  and  torsion,  the  equilibrium  equations 
can  be  derived  by  employing  the  steps  outlined  in  Ref.  1.  By  introducing 
the  Airy  stress  function  as 

N = - N + F, 
xx  xx  yy 

N = F,  (1) 

yy  xx  v ' 

N = N - F, 
xy  xy  xy 

where  Nxx  is  the  applied  uniform  compression  and  is  the  applied  torsional 
stress  resultant,  the  equilibrium  and  compatibility  equations  in  terms  of 
w and  F become 

BL^w]  - K,[F]  - - L[F,  W]  + 5^  («,„■•.;„>  - 2\y  (“•xyh'V 

(2) 


- p = 0 


L,[f]  + L [w]  + \ L [w,  w+2w°]  + -jp  - 0 (3) 

d q L 

where  the  operators  Lh>  Lq  and  Ld  are  defined  in  Ref.  1,  and  p is  the  applied 
pressure. 

The  expression  for  the  total  potential  is  given  by 

«I  - XT'  I.  + + B3F-xxF-yy  + V- J dA 

xx  A 
P 


+ ? f uy,w?  + o„w^  + a0w,  w,  + a.w?  j dA 

2 J l 1 yy  2 ’xx  3 ’xx  ’yy  4 ’xyJ 


xx  A 
P 


S [>i%v  + S3FJ  “ * S21- J»4F-xy  dA 


xx  A 
P 


- J pvdA  + (e^  + B48F  ) - *>RL  (eAVSxx  + VAV*  ) (4) 


where  the  coefficients  or ^ and  are  defined  in  Ref.  1,  and  the  symbols 

e and  v denote  the  average  end  shortening  and  average  shear  strain 
AV  AV 

respectively.  The  mathematical  expressions  for  these  quantities  are  given 


2ttR  L 


■ ¥«  - [alF’yy  + a2F,xx  + a3w> 


0 0 


+ a4W’yy  " 2 W’x  (W’x  + 2w<0]  dxdy 


2N  . 2ttR  L r 2F,  o -j 

v = -22 — 3; — f r — — y — + w,  w,  +w,  w,  +w,  w,  Idxdy 

YA V (l-v)E_  2ttRLJ0  Jq  L(1-v)Exx  ’x  ’y  y x x yJ 


Similarly,  the  expressions  for  the  end  shortening  and  shear  strain 
at  y * 0 are  given  by 

e(y-°)  = alSxx  ' L IQ  [alF’yy  + a2P'xx  + a3W’xx  + V’yy 
" 2 W’x  (W’x  + 2W0J  dX 


v(y=°)  = (T^v)e 


- L r 2F, 

: If  r SL_ 

L j _ L(1-V)E 
xx  O ' xx 

P P 


o 0~i  , 

+ w,  w,  + w,  w,  + w,  w,  dx 

’X  »y  ’x  y X yJ^Q 
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The  boundary  conditions  are  developed  in  a manner  similar  to  Ref.  1, 
and  the  general  computer  program  is  written  so  as  to  accomodate  any  combina- 
tion of  transverse  and  in-plane  boundary  conditions  (SS^,  CC^,  FF^,  i = 1,2, 3, 4) 


SS  w = M =0 

xx 

CC  w = w,  =0 

x 

FF  Q*  = M =0 

X XX 


1.  F,  = F,  =0 

xy  yy 

2.  F,  = 0;  u = C 

xy 

3.  v = C;  F,  =0 

yy 

4.  v = C;  u = C 


where  C is  a constant,  and  the  conditions  in  u and  v may  be  expressed  in 
terms  of  w and  F as  in  Ref.  9. 

The  first  step  in  the  methodology,  employed  herein  as  well  as  in  Ref.  1, 

is  to  reduce  the  governing  equations,  Eqs.(2)  and  (3)from  a system  of  coupled 

nonlinear  partial  differential  equations  to  a system  of  coupled  nonlinear 

ordinary  differential  equations.  This  is  accomplished  by  employing  the 

following  separated  form  for  w and  F. 

K 

w(x,y)  = ^ [a^x)  cos  + B^x)  sin 

i=0 

(8) 

2K 

F(x,y)  = ^ [^(x)  cos  + D1(x)  sin 

i=0 

In  addition,  if  one  considers  the  imperfection  to  be,  in  general,  asymmetric 
and  the  applied  pressure  position  dependent,  then  similar  expressions  may  be 
employed  for  w°(x,y)  and  p(x,y). 
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«°(x,y)  = Y [Ai^x)  cos  + B°(x)  sin  ^r51] 


i=0 

K 


(9) 


p(x.y)  = I [p[(*)  cos  ip  + Pi(x>  sin 

i=0 

f 

The  reduction  to  ordinary  differential  equations  is  accomplished 
through  the  following  steps. 

(1)  First,  Eqs . (8)  and  (9)  are  substituted  into  the  compatibility 
equation,  Eq.  (3).  Then,  by  employing  trigonometric  identities  involving 
products  and  the  linear  independence  of  the  sine  and  cosine  terms,  (4K  +1) 
coupled  nonlinear  ordinary  differential  equations  are  obtained.  This  sub- 
stitution clearly  shows  why  the  summation  in  the  F(x,y)  expression  is  from 
zero  to  2K.  These  equations,  which  relate  the  A^x) , B^x),  C^(x) , etc.,  are: 


for  i = 0 


11 


< • tt.  Ka  - f + ~2  l >2  [(*i + 2aj)  Aj + (bj  + 2bj)  BjJ  <10> 

R j=l 


for  i = 1 , 2 , 3 , . . . 2K 
( a ) from  the  cosine  terms 

2 


dnc7  - (¥)  < + «22  (¥)  c1  + e1[ql/1"  - 2q12  (f)  a'; 

+ <22  (f)4  AJ  + 6iAi/R  - T (t)2(ai  + 2A?X 


( 


V2°r)  Z {[(i+J)  6i+j  (Ai+j  + 2Ai+j)  + (2  ' ^j-i)^'^  6 1 i-j|: 
j=l 
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(A|  i- j | + ^1  i- j | )JAj  + 6i+j  (Bi+j  + 2Bi+j) 


[ 


- w1-”  Vii  + 2bu-ji)J  B'j + C(i+J)6i+J  (a1+j  + <j) 


- U-jI  6|i-j|  (A|i-j| + 2A| i-j| )]J  Aj 


+ {(^)6i+3  fo  + 2b»;.)  + (2 . nf_l)ii-3i»|1-Jr 


+ 2B|i-j|)J  JBj  + [6i+j  (Ai+j  + 2A°+i) 


+ (2  - 6u-j|  (au-ji  + 2Afi-j|)]j2Aj + [‘i+3  (Bi+J  + 2B°+j) 


- "i./u-ji  (Bi'i-j| + 2Bfi-ji)J  ^ - 0 


(3)  from  the  sine  terms 


dll°'i'  ' 2d12  (¥)2  D'i  + d22(T)4  B1  + 6i  klB"i'  - 2<112  (t)  B'i 


’22  (f^J  + SiB’i/R  ' ^(f)Z(Bi  + 2Bi)  Ao 


- (ik)2  1 i<i+J)26i+j  (Bi+j  + nls) + T|i-j<i'J)26|i-il: 


(Bll-Jl + 2Bli-jl)]  Aj + ['<1+J)2f,l+J  (Al+j + 2Ai+j) 
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+ (2  - 1'f-i)<1 *-J>26|i-j|  (A|1-J|  + I )J  Bj' 

-2  [-(i+J)»l+j  (»;+j  + 2 B^)  + (2  - X 

(Bi'i-ji + 2bi1-3|)]  - {(i+j)6i+j  K+J  + < j) 

+ (Aii-j| + ^ii-ji)]  jBj + [6i+j  (Bi+j 


+ 2B 


i+j/ 


+ ^i-j6  | i-j|  (B|  i-j | + 2B|i-j|)J  J AJ  + ['  6l+j  (Al+j  + 2Al+j) 


+ (2  - 1j-l)  6|j-i|  (A|l-J|  + 2Afi-j|)]  j2bj  ’ ° (12) 


where  the  prime  denotes  differentiation  with  respect  to  x and 


0 X > K 

6i  = I ; 

1 1 1 X s K 

(2)  Second,  Eqs.  (8)  and  (9)  are  substituted  into  the  equilibrium 

equation,  Eq.  (2),  and  the  error  is  made  orthogonal  to  cos  and  sin 

for  i = 0,1,2...K.  This  is  a Galerkin-type  procedure  with  respect  to  the 

circumferential  direction.  The  vanishing  of  the  (2k+l)  Galerkin  integrals 

leads  to  the  following  system  of  (2k+l)  nonlinear  ordinary  differential 
equations  (equilibrium) 


-1  X < 0 

T|x  = { 0 X = 0 (13) 

1 X > 0 
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for  1 = 0 


Ki  + 4«n)  + < (fe)  + AC  (x-)  ♦ (*C  ♦ 0 
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R d 


11. 


Po  • 


(is)2 


j=l 


+ KAj + 2aj)  Aj + (bj  + 2bj)  »)’ + (Bj' + 2bj")bj  + + 2b!'K] 


+ Kaj  + 2aj)  aj  + (bj  + 2B?)BJ  ■ 2 Kaj  + AP  cj  + (bj  + bj)  dj 


+ <a’  + A°')  c’  + <b'  + B°')  d'  + (a1'  + A°")  Cj  f (.;  + B»")d3]}  - 0 
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for  i = 1,2,. . ,K 


(a)  when  the  weighting  function  Is  cos 
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Note  that  the  equations  governing  the  response  of  the  imperfect  configuration 

1 2 

to  any  level  of  the  applied  loading  (N^.N^,  p^  and  p^)  for  a given  shape 
and  magnitude  of  the  imperfection  (A°,B°)  are  reduced  to  a systemof  (6k+2)  equa- 
tions, Eqs.  (10) , (11) , (12) , (14) , (15)  and  (16),  in  (6k+2)  unknowns,  Ai(i=0,l. , .K)  , 
B^i-l, 2. . .k) , <^(1=0,1,... 2k),  and  D1(i=l,2. . .k) . Note  that  through 
Eq.  (1),  has  been  eliminated  from  the  remaining  equations.  In  these  governing 
equations,  one  more  undetermined  parameter  is  present,  the  wave  number,  n. 

This  number  is  established  by  requiring  the  total  potential  to  be  a minimum 
at  the  critical  condition  (see  Refs.  1 and  9).  Because  of  this,  the  expression 
for  the  total  potential  i needed,  which  is  given  below 
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Note  that,  in  the  total  potential  expression,  Eq.  (17),  one  must  use 

Eqs.  (18)  and  (19)  for  e^v  and  Y^y*  but  Eqs.  (20)  and  (21)  are  listed  here 

because  in  plotting  load  versus  some  characteristic  displacement  or  rotation, 

they  are  used  for  these  characteristic  parameters.  The  reason  for  their  use 

instead  of  that  of  the  average  values,  Eqs.  (18)  and  (19),  is  that  better 

and  more  distinguishable  plots  are  thus  possible. 

Finally,  the  appropriate  boundary  conditions  are  also  expressed  in 

terms  of  A , , B . , C.  and  D . . 

i i i i 


3.  SOLUTION  PROCEDURE 

The  solution  procedure  employed  herein  is  a modification  of  the  proce- 
dure described  in  Ref.  1.  A generalization  of  Newton's  method,  applicable 
to  differential  equations,  serves  to  reduce  the  nonlinear  field  equations, 

Eqs.  (11),  (12),  (14),  (15),  and  (16),  and  the  appropriate  boundary  conditions 
to  a sequence  of  linear  systems.  In  this  method,  the  iteration  equations  are 
derived  by  assuming  that  the  solution  is  achieved  by  a small  correction  to 
an  approximate  solution  (initially  taken  as  the  linear  solution).  These 
small  corrections  are  obtained  from  the  solution  of  the  linearized  (with 
respect  to  the  corrections)  differential  equations.  The  ordinary  differential 
equations  are  cast  into  the  form  of  finite  difference  equations  as  in  Ref.  1. 
The  unknown  vector  Z contains  (12K+2)  elements 

{ z}  — { Aq  , A^ ....  A^  t B ^ . . . . , Cx ....  y ....  , 

it  u n ii  it  it  «<  «i  . 

Ao,Al’ ^£’*1 VC1 C2K’D1 D2K^ 

Note  that  the  second  derivatives  (of  and  Dj)  are  considered  as  inde- 

pendent elements  cf  the  vector  { z)  - 

By  using  one  fictitious  point  on  each  exterior  side  of  the  cylinder  ends, 
one  can  write  a system  of  (12K+ 2) (NP+2)  difference  equations,  where  NP  denotes 
the  number  of  mesh  points.  This  system  of  difference  equations  can  be  solved 
by  the  special  algorithm  reported  in  Ref.  10,  when  a unique  solution  exists 
for  a given  set  of  applied  loads.  When  the  set  of  applied  loads  corresponds 
to  a critical  condition  (limit  point),  a unique  solution  does  not  exist  and 
thus  the  solution  of  the  system  of  difference  equations  falls  to  converge. 


On  the  basis  of  these  observations,  the  description  of  the  solution  procedure 
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to  this  imperfection,  and  compare  it  to  the  critical  curve  that  corresponds 

to  perfect  geometry  and  is  obtained  through  linear  theory).  Through  the 

use  of  linear  theory  the  critical  curve  can  be  found  for  a given  structural 

configuration  (see  Fig.  1).  The  intercepts  denote  critical  loads  under 

individual  application  (from  the  typical  curve  of  Fig.  1 these  values  are: 

N = 5A  and  N = 8A  where  A is  some  numerical  constant) . From  the 

xx  xy 

cr  crL 

nonlinear  procedure,  outlined  herein  for  individual  load  application  and  a 


given  imperfection,  one  can  find  the  nonlinear  theory  critical  loads  (see 


Fig.  1,  N = 3A,  and  N 

xx  xy 

CrNL  crNL 


= 7A) . In  order  to  find  the  complete  curve 


one  can  (a)  fix  N = A and  2A  and  employ  the  methodology  by  increasing 
xy 

N to  find  points  I and  II  and  (b)  fix  N = A and  increase  N to  find 
xx  r xx  xy 

point  III. 


There  is  an  alternate  approach  in  constructing  the  nonlinear  critical 
curve.  This  approach  requires  that  both  Nxx  and  N be  increased  linearly 
(along  line  OP  corresponding  to  some  angle  9)  in  order  to  find  the  critical 
combination,  pointlll,  for  that  9 value.  Then  vary  9 from  zero  to  rr/2  and 
the  complete  curve  is  generated. 

Both  approaches  are  incorporated  into  the  computer  program  (Appendix  A) . 

Regardless  of  the  approach  used,  the  ratio  ( 01 I I | / [ OP | is  a measure  of 
the  imperfection  sensitivity  (instead  of  j QIII  | / |QP'|).  This  ratio  may  be 
called  knockdown  factor  and  denoted  by  A. 


I 
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Fig.  1.  Critical  Conditions  Under  Combined 


4.  NUMERICAL  RESULTS  AND  DISCUSSION 


The  present  methodology  is  demonstrated  through  a number  of  illustrative 
examples.  Numerical  solutions  are  obtained  by  employing  the  Georgia  Tech 
high  speed  digital  computer  CDC-CYBER  70,  Model  74-28. 

A general  computer  program  is  written,  which  includes  the  following 
desirable  features  (Appendix  A): 

(1)  It  is  applicable  to  a stiffened  configuration,  in  either  or  both 
directions,  as  well  as  to  an  unstiffened  configurations. 

(2)  It  accomodates  all  possible  boundary  conditions  (SSi,  CCi,  FFi,  etc.) 
and  it  can  easily  be  modified  to  accomodate  elastic  end  restraints. 

(3)  The  number  of  Fourier  terms,  K,  can  be  as  large  as  needed  for 
accuracy.  The  same  is  true  for  the  number  of  mesh  points,  NP,  in  the  finite 
difference  scheme. 

(4)  The  shape  of  the  geometric  imperfection  is  unrestricted. 

(5)  It  is  applicable  to  any  individual  or  combined  application  of  uni- 
form axial  compression,  torsion,-  and  space-dependent  lateral  pressure. 

(6)  The  CPU  time  required  to  obtain  a solution  is  reasonably  small.  For 
example,  by  using  K = 1 and  NP  =*  57  (826  unknowns)  a solution  (critical  load 
and  all  intermediate  steps)  is  obtained  in  18  seconds.  For  low  load  levels 

a convergent  solution  is  obtained  through  two  iterations.  For  load  levels, 

approaching  the  limit  point,  a convergent  solution  is  obtained  through  six 

iterations.  The  solution  has  converged  if  *-he  percent  difference  in  response 

-4 

between  two  consecutive  iterations  is  smaller  than  10 

The  numerical  results  for  all  illustrative  examples  (ten)  are  presented 
in  tabular  form  in  Table  1.  The  geometric  parameters  and  load  conditions 
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for  each  example  number  are  shown  on  this  table.  In  addition,  the  perfect 
geometry  (linear  theory)  results,  the  imperfection  analysis  results  (present 
nonlinear  theory),  and  the  associated  knockdown  factor.  A,  are  shown.  For 
all  ten  examples  the  boundary  conditions  are  taken  to  be  classical  simply 
supported  (SS3) . Poisson's  ratio  for  examples  4 and  5 is  0.3333  while  for 
the  remaining  examples  is  0.3. 

Examples  1,  2,  and  3 have  been  reported  in  Refs.  1 and  7.  The  imper- 
fection considered,  for  these  examples,  herein  is 


°f  \ _ mrrx  ( 

w (x,y)  = t sm  — 


ny\ 

R / ‘ 


Although  a more  general  imperfection  shape  is  used  in  this  paper  (in  Refs. 

1 and  7 the  imperfection  shape  was  taken  to  be  symmetric,  cos  ^ only)  the 
results  are  the  same  because,  in  the  absence  of  torsion,  there  is  no  coupling 
between  sine  and  cosine  terms. 

Examples  4 and  5 correspond  to  an  unstiffened  configuration  under 
torsion  only,  and  they  correspond  to  cases  that  have  been  worked  out  previously 
(see  Ref.  5)  by  employing  a Koiter-type  analysis.  The  difference  between 

4 

these  examples  is  the  curvature  parameter  (Z  = 10  for  Example  4,  and  Z = 10 
for  Example  5).  A comparison  between  the  present  results  and  those  of  Ref. 

6 shows  very  good  agreement  for  both  cases.  The  imperfection  for  these  two 
examples  was  taken  to  be 

5 

w°(x,y)  = 6 £ [Am  cos  ^ + Bm  sin  sin  ~ (22) 

m=l 

with  6 = t and  where  A and  B are  the  elements  of  the  corresponding  perfect 

m m 

geometry  linear  theory  eigenvector. 
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Note  that  in  the  case  of  torsion  there  are  two  different  eigenvectors 
(linear  theory),  one  corresponding  to  positive  torsion  and  one  to  negative 
torsion.  Since  the  imperfection  shape  is  considered  to  be  similar  to  this 
eigenvector,  the  analysis  is  performed  to  check  the  effect  of  one  on  the 
other.  It  is  found  that  an  imperfection  shape  similar  to  the  positive  tor- 
sion eigenvector  has,  virtually,  no  effect  on  the  load  carrying  capacity 
when  the  torsion  is  applied  in  the  negative  direction.  (The  shell  is  in- 
sensitive and  thus  A = 1) . This  reinforces  the  contention  that  the  shell 
is  sensitive  to  imperfection  shapes  which  are  similar  to  the  corresponding 
perfect  geometry  linear  theory  eigenvector.  This  observation  is  made  for 
all  other  examples  with  torsional  loads. 

Examples  6 through  9 correspond  to  a stiffened  configuration  (same  for 
all  four),  imperfection  shapes  characterized  by  Eq.  (22),  and  different 
imperfection  amplitudes  (6  = t,  2t , 3t,  and  5t  respectively).  The  results 
for  these  cases  are  also  presented  graphically  on  Figs.  2 and  3.  Two  ob- 
servations are  worth  mentioning  here.  First,  the  "guessed"  postbuckling 
behavior  (see  dotted  line  on  Fig.  2)  indicates  that  shells  loaded  in  torsion 
are  not  as  sensitive  to  geometric  imperfections  as  those  loaded  in  axial 
compression.  This  Is  in  agreement  with  the  qualitative  results  of  Ref.  3. 

Second,  stiffened  configurations  of  the  same  curvature  parameter  (based  on 
actual  thickness)  or  a weighted  curvature  parameter  (based  on  a weighted 
thickness  that  includes  the  smeared  stiffener  contribution)  are  not  as  sensi- 
tive as  the  corresponding  unstiffened  configurations.  According  to  Ref.  6 
(see  Table  1;  Case  I,  and  Fig.  6 of  this  reference)  the  knockdown  factor 
for  an  unstiffened  configuration  with  Z < 96  is  definitely  smaller  than  0.6 
(A<0.6),  while  for  these  stiffened  geometries  AS  0.86. 
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Finally,  Example  10  corresponds  to  the  same  stiffened  configuration 


as  Examples  1 and  6 through  9 but  under  combined  application  of  axial 
compression  and  torsion.  The  imperfection  shape,  for  this  example  also,  is 
taken  to  be  similar  to  the  classical  perfect  geometry  eigenvector.  Since, 
for  this  geometry,  the  analysis  is  performed  for  individual  load  application 
(Examples  1 and  6) , the  knockdown  factor  is  computed  as  described  in  the 
previous  section  (see  Fig.  1). 


0 4 8 12  16  20 

y (y=  oi  x i o2 

♦ 

Fig.  2.  Load  Versus  Angle  of  Twist  for 
Imperfect  Stiffened  Cylinders 
(Examples  6-9). 
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Fig.  3.  Knockdovm  Factor  Versus  Imperfection 
Amplitude  (Examples  6-9). 
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5 . OPTIMAL  CYLINDERS  IN  TORSION 


Since  no  stiffened  or  unstiffened  shell  configuration  is  free  of  ini- 
tial geometric  imperfections,  the  optimization  of  such  configurations  re- 
quires the  specification  of  both  the  amplitude  and  shape  of  the  geometric 
imperfection  and  the  use  of  a nonlinear  buckling  analysis  in  the  solution 
methodology.  Both  requirements  present  difficulties  of  various  degrees,  the 
most  serious  of  which  is  the  inclusion  of  the  nonlinear  buckling  analysis  in 
the  optimization  algorithm.  Because  of  this,  an  alternate  approach  will  be 
employed.  This  approach  is  much  simpler  in  terms  of  execution  computer  time 
and  does  lead  to  a reasonable  final  solution.  Before  presenting  and  demon- 
strating this  alternate  approach, let  us  state  exactly  the  problem  in  hand: 
Given  an  internally  stiffened,  thin,  circular,  cylindrical,  imperfect  shell 
of  specified  material,  radius,  length,  and  imperfection  (amplitude  and  shape) 
find  the  size,  shape  and  spacings  of  the  stiffeners,  and  the  skin  thickness 
such  that  it  can  safely  carry  a prescribed  torsional  load  with  minimum 
weight. 

The  alternate  approach  consist  of  the  following  steps: 

1)  Employ  a safety  factor  with  the  applied  load,  which  accounts,  among 
other  things,  for  the  imperfection  sensitivity  of  the  stiffened  shell. 

2)  Perform  an  optimization  of  the  corresponding  perfect  geometry  con- 
figuration by  employing  the  methodology  described  and  demonstrated  in 
Ref.  11. 

3)  Using  the  specified  imperfection,  perform  a nonlinear  stability 
analysis  (see  chapters  2 through  A)  on  the  optimum  configuration  (from 
step  2).  Through  this  analysis  find  the  corresponding  knockdown  factor. 
If  it  is  close  to  the  guessed  number  used  in  step  1 proceed  with  the 
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remaining  steps,  otherwise  modify  the  safety  factor  appropriately 
and  repeat  steps  one  and  two. 

4)  On  the  basis  of  step  2 the  design  space  surrounding  the  optimum 
configuration  is  available.  Choose  design  points  in  this  space  (their 
weight  is  higher  than  that  of  the  optimum)  and  perform  a nonlinear 
stability  analysis  in  order  to  establish  the  corresponding  knockdown 
factors.  This  comparison  will  establish  the  optimum  design  point 
in  the  presence  of  specified  initial  geometric  imperfections. 

The  entire  procedure  is  demonstrated  herein  by  employing  a design 
configuration  of  Ref.  11. 

6 . 1 Design  Example 

The  example  used  herein  to  demonstrate  the  optimization  procedure  is 
Example  1 of  Ref.  11.  The  specified  parameters  are  given  below: 

R = 85  in.  ; L = 100  in.  ; E = E = E = 10.5  x 106psi. 

x y 

v = 0.33  ; psk  = Px  = Py  = °-101  lhs/in.3 

The  applied  torsion  (stress  resultant)  is  taken  to  be  300  lbs/in.  and  by 
assuming  a safety  factor  of  approximately  1.4  to  account  for  imperfection 
sensitivity  (knockdown  factor  = 0.715)  the  design  applied  torsion  used  in 
the  linear  buckling  analysis  optimization  procedure  (Ref.  11)  is  418.5  lbs. 
In  considering  various  types  of  stiffeners  in  Ref.  11,  it  is  concluded  that 
the  best  combination  of  shapes  corresponds  to  T-stringers  and  rectangular 
rings  (TS-RR) . 

The  imperfection  shape  assumed  for  all  configurations  is  taken  to  be 
similar  to  the  linear  theory  buckling  mode  and  the  amplitude  of  imperfec- 
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tion  is  taken  to  be  the  same  for  all  design  points  employed  in  the  com- 
parison (6  = 0.1215  in.). 


Note  that  the  imperfection  is  characterized  by 


w°(x,y)  = 62  (a  cos  ^ + B sin  sin  ^ 
« \ m R in  R / L 

tn— u 


where  Am  and  Bm  are  elements  characterizing  the  linear  theory  buckling 
mode,  each  divided  by  the  coefficient  of  the  most  influencing  term.  Through 
this  the  largest  of  A^,  is  equal  to  1.0  and  the  remaining  absolute  va- 
lues are  less  than  1.0.  Thus,  the  maximum  amplitude  of  the  imperfection 
is  given  by 


6 = 6 


I i {* 

m=0  ' 
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For  all  configuration  considered,  6 is  varied  such  that  6 is  always  equal 
toO.1215  in.  This  represents  an  amplitude  of  approximately  2.5  times  the 
thickness  of  the  linear  theory  optimum  point  (t  = 0.05  in.).  The  under- 
lying thought  here  is  that,  in  order  for  the  comparison  to  be  meaningful, 
the  amplitude  of  imperfection  must  be  the  same  for  all  configurations. 

The  geometry  of  all  configurations  employed  in  the  comparison  are 
given  on  Table  2.  The  symbols  used  are  the  same  as  those  of  Refs.  1,7,8 
and  11.  Point  1 corresponds  to  the  optimum  point  as  obtained  by  employing 
linear  buckling  analysis  in  the  optimization  procedure.  Points  2 and  3 
correspond  to  the  lightest  configuration  for  a given  skin  thickness  or  cur- 
vature parameter,  Z.  Points  la  through  le  correspond  to  designs  surround- 
ing the  optimum  point  but  for  the  same  value  of  the  curvature  parameter 
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as  the  optimum  point.  Points  4 and  5 are  discussed  later. 

6 .2  Discuss  ion 

All  of  the  design  points,  discussed  above  and  listed  on  Table  2, 
represent  the  linear  buckling  theory  optimum  design  (pt.  1)  and  designs 
surrounding  this  optimum  (pts  2,3,  la-le) . Each  one  of  these  points  was 
then  analyzed  in  order  to  find  its  load  carrying  capacity  (buckling  load) 
in  the  presence  of  the  same  amplitude  imperfection.  The  linear  (perfect 
geometry)  and  nonlinear  (imperfect  geometry)  buckling  loads  are  given  on 
Table  2. 

The  goal  of  the  entire  procedure  is  to  find  the  lightest  possible  con- 
figuration, which  can  carry  safely  a torsional  load  of  300  lbs/in.  A com- 
parison among  all  designs  shows  that  designs  1,  and  la  through  le,  all  can 
carry  safely  approximately  300  lbs/in.  The  lightest,  though,  of  all  these 
designs  is  point  1.  Therefore,  one  may  conclude  that  the  optimum  design 
of  the  imperfect  stiffened  cylindrical  shell  is  given  by  design  point  1. 

In  order  to  strengthen  the  above  conclusion  one  more  question  must 
be  answered.  Before  posing  this  question,  though,  let  us  state  two  obser- 
vations based  on  linear  theory  optimization  procedures.  Given  a stiffened 
cylinder  of  specified  radius,  length  and  material,  first  by  increasing  the 
value  of  the  applied  load  (small  increments)  the  corresponding  optimal 
weight  increases;  and  second,  the  same  small  increments  in  applied  load 

result  in  small  variations  in  the  optimal  geometry  (a  ,a  ,\  ,p  ,p  , 

x y xx  yy  xx  yy 

and  Z).  Having  made  these  observations  (true  for  other  load  cases  as  well) 
the  question  that  must  be  answered  is  the  following:  Is  it  possible  to 
decrease  slightly  the  applied  load  (by  increasing  the  knockdown  factor) 
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and  have  an  optimal  configuration  which  is  less  imperfection  sensitive,  thus 


resulting  in  a lighter  configuration,  which  can  carry  safely  N = 300  lbs/in.? 

xy 

Also,  how  imperfection  sensitive  is  an  optimal  configuration  obtained  from 

the  linear  buckling  analysis  optimization  with  a slightly  higjier  applied 

load?  In  order  to  answer  these  two  questions,  optimal  configurations  are 

generated  for  N = 394.7  lbs/in.  and  435.0  lbs/in.,  and  presented  as 
Xycl 

points  4 and  5 on  Table  1.  As  expected  the  corresponding  weights  are 
slightly  lower  and  higher  than  that  of  point  1,  respectively.  By  performing 
a nonlinear  buckling  analysis  on  these  configurations  (pts . 4 and  5)  it  is 
clearly  seen  that  both  are  more  imperfection  sensitive  than  that  of  point 
1 and  thus  we  may  conclusively  state  that,  the  optimum  geometry  for  the 
problem  in  hand  is  given  by  design  point  1.  Note  that  the  last  row  of 
Table  2 gives  the  ratio  of  total  weight  to  the  critical  torsion  (shear 
resultant)  for  each  design.  This  ratio  can  be  thought  of  as  a performance 


index  for  each  design;  the  lower  the  ratio  the  better  the  design. 


TABLE  2.  Performance  of  Various  Design  Configurations 
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COMPUTER  PROGRAM 


BLOCK  DrAGRAhA  . 


Si^ce  4U  IS  CXva  ovy^d  VwO  JL'^W 

"4Us.  <rw£  e_w\ y^Io^c-S.  kv\  {.7.  > v^ouJbe/V  is  ve | e*t~>re ' 

4o  ife/yi  Cw\A  ^ppev\JU*  S o^_  ^”2.* 


fV\A\W  PfcOGrR6fY\ 


Dfti-os 


xmpEC.yr 

XvApe^r^etlC«»i 
<Vv\vj  poiwt^l^v 
(Xmvj  PovAY^er  iervv^ 

'bVA.Iiro  i\inv\€. 

W'fitiev\  V)Y)4l'« 

user. 


POTTERS 

A^gorrVKvw  'V0^ 
i>oQ.u»vi^  4i%e  in 
y>  eotv  ^ vj  srt&wi 


CcxicvAi^es  4iv€  i^wui-  ipov^L  QxmxJL  cavA 
-Vta,  WcHiw^  vno^e  •,  pt'iJu  o Jt  oJLi 

(d^A.  <X/v4  OVA^piA^  pOsASKmelxvs 


Colcuiules  (jJi 

Coe^ic\ev\Vs 
ujLacU  okvc  t»»<le- 
^ev\le>A, 

C»  vCvavw 

UJft.u£  YMHw\la£-Y' 
Cv\), 


C0 E^NN 
Cq1ua<UJCa>  saW 

oo  \Ca  caJu 

LoUicJU  OYY«  dcy«»4 

€vrt  OVA  +L«- 
Ca  'TCAav*  ftl 

UUtt-AJC  Y\UAv\iev 

O}. 


TC-f^F 
CreeJti-,  nectars 


.Bm  - 


-c*'d" 

w\ 


I 


r > 

4-6 

v y 


y > 

v4y 


POTSN 

G^CvJLOut-,  4io.  pci' 
ewiicj  eAAC.N'^  ,€mJ1 
^IvxrriiiMtvx,^  <fc/wJ[ 
$Vv^v\ 

€ay>W,«(v°);YM 


o^-pireuw>vAi  sotu- 
ti«v\  w>  O'fi^^or'io 
© ii<XA/v\  v\«  4-  tfe-j 

Y'ftfcid-n  (xvi Jt  - 

YVAOY-t^es  ti  t)^ 

GOvVxvW  0 W *»fc»sijL- 

^e^X. 


S eJts  ALe  ^aiuA 
^ iUe  pcvAv  caa- 
Gwv  XocdL  ocktx(, 
^<,€€  ^£*VJ  civ<wi  I 
yrspftmd^ 
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fXfeCGi 

ces  (Xc«-0'r<X>  v\c^  "Go 
Pottexs 


r > 

42 

v y 


vsvivw 

Modv<x  ofe*i»tv«n 
A=fc-C*D 
C^ee  &4-  4^) 


r \ 

4-5 

v y 


yn\y 


(V\o^,«-tx  <Jjp 

A -B*C 
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v y 


INVERT 


1 G>0«S  'V&'C 
ww-erse  ©V 

t 

i 4U?  vviodtri  x, 


XWRITE 


Wrifce  iwW- 
Vv\e«Uak.  soULwn 
o*  divett  <xctett 

iota.  <*e4.  (21) 
fcvdl  . 


Bounds 

w\cx4v'xc 

^4vje  Oy^|irO^Tvo4<.  koVA-W^fi.- 

CcvN^4»CY\i  js'i'z.ytf  1^=§] 
{.  ^ee  ?*■£. 


XRCAT^ 


£eodu  4-l~e 
\ \n4e^wx  e JU^Jbe. 
Soft>j3C\X>v\ 

(21) 

&/\Jl  (^G)  , 


R^T(> 

Gcew^rro^jeo  wcO[\v\Cfeo 
')A  I 0 ^ Qo-wvv|<jr1a 

&/WT-  ^uiiU  Wiix.wv  ex^vofctaw 

r*]7.M<i2',Wt.j(«4-R).' 


AV_£ 


CtnJxAAXsJ^i  cocJ|^»  c\€^io  \du) 

Gc  T«uv\€«r  W.%Co  -to  . CG^.  I*), 


34 


CCM(Y\OK 

C^MrA^N/cx  mtg/N£Q> ?0T/)  (ym  Cso°V 

NtQPtfiT  — NuamW  'poivko  */*  O^uoJ.  <Uv«cturrt  . 

tV%X  (S"oo)  - Tt\ft.  crrAav  j (VTC  (4}  , <r|.Ec^,  ^ °V  ^V^-evA  ©jb^ovA-^v*. 

■2)  C0<VWV\{*n/e>O\jND/  LSI  , \-SN  , 

•M'ive-rt  O-j.  V)OUV% Uyv\^Al-^CnA^  ftX  yoiv^- 

cjl  v^2  JU/A  \po\A  ^usN)  o-\-  S-Lft  A%eW. . 

3)  C <}frMlY\^  n/ \^I  ft  Ffc/ DELTA,  ftLl  ;GA  1,  AL2,  ^Te^GAZ 

Co^tieAs  «r^  vN»\je.  •|<w  vn  ( fce^AZ). 

4)  c^A\A\^(vjy/f^ug  l/wiF^uz , ki f we, . . - . vs 

C^iY\m^rv/^ugz/\clO;k:ll; ...  fcls^Vikz. 

Fou'neA'  JUwuJt  (k=kF^U£)  (maA  ^o^ovAmAe^i 

Ae^v\jLeA/A  <rv\  V. . 

5 ) C0H\fl\0rt/<xStfiA/(llJt>l);  D14 ,0 ,2£?; 

^bUeid  ^exm^t^c  ^©^rcMAAcjUvj  \z}. 

/FACTOR/ Cl  ,C2;  - . . C 12. 

Coft||iCA€AAo  crv\  C. \ 'rOoww. -^a'&AaaJ?  oJftA/e  vavaavJI^v  £v\}, 

7)  C<^(V\iw^n/fact2/  Oli,  -^>l4  ^At,--  t)A4,DBt,-  DB4,XNi; £xxy> 
a^'oe^  , 'K1-^)ft1-a4,Voa>3j^4-,^^xxf  (*«  «4- i*). 
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*&)  c$<V\iv\(^n/fact3>  /xLjXH  . 


X L - \_ew.^tU  ^ Xw  -M* uUes\ 

S)  c^(v\(v\^k//cotsx/  xzi.(s~oi^)  (roi) 

Dw^cF.  access  ^©Jia.  seio  Q 1 O'^  '2.7, 

10)  c0(v\fY\<^N  / Pees  l /Avm  (100,3) Awmp (100,3),  aw(v\pp  (100,3). 


VecWs  , My*  ,Aw\  °V  S o^vA^jw 

11)  c0*v\rv\0N  / Peesz/^ww  (400/2),  fcwwP  (4oo,7),b>ww\pK*-o0/0. 

Vectors  Byy\  e^.  ^Ye^voas  s oQ.iAc&'v'x  . 

1Z)  C<j/>U\W\j6u /p«es  3 / CFM  (10Q;4-),CFN\?  ( 1 00^4 ),CF  (WPP  (l00,4) . 


UecUs  CW>C^  ,C^  4 ^YfcAMCVS  % O^LAoTvy 

\ 3)  C<6t\MW^N/pC.cSA/DV^(\00^);P^  (loo/O^P*^  (w°/0-  1 

vecJUvs  Oy* /Dr*  >tC  yrtAJiovAS  So'IaAa^. 

14)  C^(V\m^(y/?fiesb/AWZ(|0OJ3)>,^WZPj(lCO;3)>)AV\J%PV3  (lOD^). 

S^Vv^-tWi  c VW\^1 .Pyv^ecX^Vk  W\ojle  , l\°  ' , A°  /f. 

!§■)  c^m^w/PftesC/BWz,  (ioo,?),BwZP(ioo,£),SwzPP(ioo/0. 

^Vy\V\S»^yA\vv\ety  i c \ vv\(Ji^Ar^C>tjU>erv\  w\©  Ae  2)°;Bo/;30,/. 

1&)  C^»A<v\y)hj/^es6N4/  ftwiyfe  (?)3),ftW>w?>Pfe 

V/SectoYS  (A  ,B)  <0 1 J'revJv  OWS  So^aa1/VAt-v\  ocK  ^'cA-VlVows  Y°'Ao 

17)  (2,4),CFfVN?PB  (2)4  ),t)F(V\8  (^^O^WpPBfc/O. 

V^ctejv-S  (c.jO)  <r(  (JViaoiovAi  soQaaI^tv)  a^.  po\‘Ao. 
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TL^Af)  rT.  fov  o^Kj\^  <m  cojxJ^,  ^ 

Xl^ At)  - 3 A$si<^  W*x  If  <^>  ^®m***©*k.  cab4*fco 

Xl-^AT)  r 4 fW\'<Y*\  <X/>  ^e*vv*A/v'«v&  , WXv  <m~4  v °^e 

C«vvv^eJtxX  Ip^  ^ uowAjOwifc  (tJxx  = m*f) 

OKSlxoJL&i-^o  Hxx  ^ W'A  . 


i,ro-w> 


i?ftesv*ae  - i_ 
X ?A£S  U9e - Z 


OrwiLsw^  ^re^^vf. 
x-Ae^evX&wL  pressure 


6E  1Z.4 


cov\s"U^£  fYfisvAre 


( p (x)  ^rd.WEQP^)  6 El?  ,4 

X-^t|>eAA^vd-  |pres^i  Fav  eJU. 
^OVxaA^  \ X -Ai'rt^UvN^  J r 4 ;N6q?^T. 

Mote  , X '>f  SSu.M~e  c&A\v\/si  $«.  x -JU^«kl 

&<A.  ■j-w  ll^AO  : 4, 


jL0Afc-=\_^ i A7 
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© 

GE12.4] 

GG12.4 

X^XXjDNXX^NXY.KecaS  ,RT 

+ 

PFAC  &S>V  AC  jXNXX^KNXY,  ACCUd,  ex 

X^x>  - I y.ii»<A  0A»«^  C£ry^j)ire^Ov\ 
®NXX~1v\u»re>v\e^4  a^O^voJl  omyf«*i*» 
XNXy-^\xe^  “tbrslov^ 
ftccuft.  - ^CCutTXC^  (/O^IaS. 

^evce/vd-- 

1 9T  *V\U/)uvW>m  VmaaJUv 

v*  *4ti  Or^CASt  !lo<xA  O^tAA- 

i , 

t 

?£AC  " futlrT  4^  jyyfSSAATl 

"t/A  y^A-ioJl  o^aJL>tXa-*vvo 

DPFAC'  Ivn^w^A-  w>  <&«*  ^a.ii*r 
XNXX  — (Xi>\o5  Cwyv^'|>Tfev;vffv\ 

VNXy  ~\»TScev\ 

Acc  uft  J C-x  — “sowt  &o  »*  (7S) 

I 

1 

XKxy^w^v,XNxxj(\ccuft1  ex 

XNxy  - Iw.UoJ  ta>vs\  <jv\ 
t)N^y  - i.y^ore>v\evt  ^'Wmoyn 
XNX>  - V^XeA  Oow^r€SSi(rv\ 

^CCU(? 
fci 


0* i | v\  © 


6 €12.4 

PFAC,  DPfrPyC  ,?WX  ,VNXVf  Acc  ua,fci 

9 FAC  ~ V-a. Jbrr  V'AvJU^yi^  ycesswre 

■*«  CXmJUVitfv'i 

\)?FAC  - ^wrt wei  u (A<rvf  ^dUc 
9ni  - V®vii«v  uWd*  vwJLL^JUto 

Acc^fcT  U (ft?) 
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NNM,  LNNvi  ;1LNW  ) 

t'lNN  (v>)  TU  ^ro^vwwv  4it  (U*»t  ^>r4U^  y^, 

LNNU  - 0 tL?  jJV^K^YOjWv  Ao4rT>  v\xA/  civ€xV-  v>)lv'Jlx  V\  -V>  oSaa-6» 
Qjpv^O  yvvivNxvv\w.  ^Atv^tv*J)  ( toWi^  ) 

LnNN  - 1 'TU  ^AAAaI^O  A^ajI;  i/\*wb  yo\vJu  ^*Y  NilWy 

<j>vijuJUJieA  N-Ab<£  jutaJi  ^Ajc^Vv«J)  -^w  V\-  ^AaJbv*,-) 

Vvory eyw^U w(^  VJ(AN  (Xa*4^  Q/£«vx)l  j^p\kH/Y 

v-Vbi  Q/yw*A  \ vOc  ^ Vi  N N r 

LNNN  r Z.  ^TW  jiyoC^v'O^'  • CxjSLu*.^1«JU^  4ist  jfcboJL 

^ &\&v\biC>J]  <Sv\a^  oA  'U/4  \V\CWJ1  ^^vv'VU* 

T \^  (V] W i,  W|^6A  JLwV\b  <JV\  V\  -^W  yj  Vva  (A\ 

^)Y«^vWv  coJIcaaJU^/3  AasWfi  ^oWA^a*, 

\AA  (rvlfi^  bo  |>WvJl  uAvCtby  WY  ^p^<2A 

Wv  \ v\  \ W V/VWV.  ^ oV'i/V'b^aX.  , 


U C <25  Nf~T  N 
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IC^NTM  ^ i.  TC^y*  Vo  GYvvcAi^ 

"to  ^‘VA^V\1 

UC^WTM  ^ A.  'tCwj  \o  4b<  Ibvvb 


PROGS AM  MA IN  (INPUT, OUTPUT , T APE5  = I NPU T , TAPE 6=0UI PUT. TAP E21» F APE22) 
C . MAIN 

C POST  BUCKLING  OF  STIFFENED  CYLINDRICAL  SHELLS  UNDER  COMdINEO  LOAD  OF 
C UNIFORM  AXIAL  LOAD,  LATcRAL  LOAD  AND  TORSION 

• COMMON/FOUR1/KFOUR,K1, K2, K3,  K4 , K5  , K6 , K7 , K8  , K9 
- COMMON/ FOUF.2/K1  0 , K 1 1 , K 12 , K 1 3 , K 1 4 , K15  , KK2 
— CUMMON/CINTG/NEQPOT.MI (501) 

- COMMUN/COIS</I21 (501) ,122(501) 

- COMMON/FIOFR/OELTA,ALl ,GA1 ,AL2,OT2,GA2 
COMMON/PRESl/ANM (100, 3) . AKMP (100, 3) , AWMPP( 100,3) 

COMMON/PRE S2/8WM  (100,2)  , dWNP  ( lu  0, 2)  , dWMPP  (100,2) 

CCMMON/PRES3/CFM (100,4) ,CFMP ( 1 0 0 , 4) , CFMPP  (10  0,4) 

COMMON/PKES4/OFM ( 1 0 0 , 4 ) , OF MP ( 1 0 J , 4) , DFMPP ( 100 , 4) . 

CUMMON/PRES5/AWZ (1C0,3),AK2P(100,3>,  AWZPP ( 10  0, 3) 

COMMON/PRE S5/OW  Z(100,2),8WZP(100,2),  8WZPP (100,2) 
C0MM0N/RES6N1/AWMBI2, 3) , AWMPP8 ( 2, 3) , DWM3( 2 ,2) , BWMPPO(2 , 2) 
COMMON/RES  GN2/CFM3  (2 , 4 ) , CFMPPB  ( 2 , 4 ) , 0Ffi3(2 ,4) , DFMPP3(2 ,4) 

- COM MO N/GEOM/RR , DD, HI 1 , HI 2 , H2 2, Q 11 , Q1 2, Q22 , Dll , D 12 , D22 

- COMMON/FACTOR/C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12 

— COMMON/FAC  T 2/OL 1 , 0L2 , CL  3 , DL4 , D A 1 , OA2 , 0 A3, D A4 , 082, Db3 ,0  84 , X NI , EXXP 
—COMMON/F  ACT  3 /XL , XH 

COMMON/XXL OA D/ A XPRES( 100 ,3) , BXPRES ( 1 00 , 2) 

- COMMON/ BOUND/L S 1 ,LSN 
DIMENSION  TI  (10)  ,xyPRES(100,3) 

DIMENSION  A? (52, 52), dP (52. 52), CP (52, 52), PR (52, 5 2) ,GP(5 2,1) 
DIMENSION  XP (52 , 1) , Tl (52) ,CC (52 ) ,MT (52) ,V1 (2704) 

C 270^=52*52 

DIMENSION  ARUM ( 3) , BWWM ( 2 ) , CCM ( 4 ) , OQM (4) 

DIMENSION  ANCON (20,3), BWCON (20,2), CCON (20, 4), DD ON (20, 4) 
EQUIVALENCE  (AP  ( 1 , 1 ) , V 1 ( 1)  ) 

CALL  OPENMStZl, 121,501 ,0) 

CALL  0PENMS(22, 122,501,0) 

ECONv/  = 0.  0000  01 

■ IIA  X N=  52  • . ' 

MAX 2= MAX  N*MAXN 
NRHS  = 1 
!(J=  10  0 
NAW  = 3 
N B W = 2 

NF=4  ' ' " 

M3=ioo 
M4  = 3 

1111  WRITE (6, 20) 

REAO(5, 10)  (Tl(I) ,1  = 1,9) 

WRITE (6,60) 

WRITE  (6,10)  (TI ( I), 1 = 1, 9) 

RE  AD ( 5 , 100)NEQPOT,KFOUR, LSI,  LSN , lPRI NT , LMOD 
IF (LPR1NT. EQ.l) LMOD=l 
I DDE  T =1 

RE AU(5, 200) RR, XL,XH,ELAS»XNI 
READ  (5, 2 GO)  XLAMO  , YLAI-ID  , E X,  EY  , RHOX  , RHOY 

Q *♦♦»***♦*♦*♦»♦»***  **+**»*, ♦****»♦******♦#< 

CALL  COEFF (EX, EY.XLAMD.YLAMD, RHOX, RHOY, ELAS) 
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WRITE  (6, 30  0) NEQP0T,KF0UR,LS1 ,LSN 
W p I T L ( 6 , 4 0 0 ) RR , XL,XH,ELAS,XNI,UD,EXXP 
WPlTt (u, 5Z31XLAM0, Y L AMO, EX , E Y , RHOX , RH0Y 

CALL  IHPERF 

WRITE  (6,5081 
J1  = 0 
WRITE (6,510) J1 
WRITE(6,520) 

X X = 0 . 

□0  85  Il=l,NEQPOT 

WRITE  (6,50 9) II,  XX,  AWZ(  11,1)  , AW ZP(  II,  1)  , AW ZPP.( 1 1 , II 
XX  = XX  + OELT  A 
35  COMTINUE 

□0  86  Jl=l,<FOUR. 

WRITE (6, 510) J1 
WRITE (6,520) 

XX=U . ' 

DO  86  11  = 1 *NEQP  OT 

WRITE (6, 5091)11 ,XX,AWZ(I1, Jl+1) , A WZP ( I 1 , J 1*1) , A WZPP ( I 1 , J1 * 1 ) , 
ldWZ(Il,Jl)  ,dWZP(Il,Jl)  ,BWZPP (II , Jl) 

XX=XX*DELTA 
86  CONTINUE 

IF (LPRINT.NE.l) GO  TO  39 
WRITE (6, 50 01  DELTA, AL1,GA1,AL2, BT2.GA2 
WRITE  (b,  5011H11 , H12.H22,  Qil,ai2,GL22 
WRITE (6,502) Dll , 012 , D22 , 062 , DB3 , DB4 
WRITE  (6,503)DL1,UL2»DL3,D14 
WR I TE (6, 50  A) DA  1 , DA 2,0 A 3, OA A 
39  CONTINUE 

DO  63  Il=l,NEQPOT 
DO  64  Jl=l »Kl 
A w :•)  ( 1 1 , j l ) = 0 . 

. AWUP( II, J1)=0. 

AWMPP  (I1,J1)=0. 

AXPRcS(Il,Jl)=0. 

XXPRESUl, J1)=0 . 

IF ( Jl .EQ.Kl) 00  TO  64 
BWrHIl,  J1>  = 0. 

DWDP( II, J1)=0. 

OWMPP ( II , Jl) =0 . 

BXPk£S(I1,J1)=0 • . 

64  CONTINUE 

DO  65  J1=1,KK2 
* CFM(I1, Jl) =3. 

CFNPC II , Jl )=0. 

CF.1PP  (II  ,J1)=0  . 

DFMCI1, J1)=0.  • 

OF  NP ( 1 1 , Jl > = 0» 

0F«PP(I1 ,J1) -0. 

65  COMTINUE 
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6 3 CONllNUc 

READ(5«100) ILOAD, IP RE SS 
IF(IPRESS.£Q.2> GO  TO  637 
R£  A D ( 5 , 2 0 0 ) Pll 
JO  639  I1=1,NEQP0T 
63 9 XXPRES  (II,  1)=P11 
GO  TO  b40 

63  7 READ  (5,203) ( XXP  RES  ( 1 1 , 1)  , I 1=1 , NEQPOT) 

115  640  CONTINUE 

WRITE (6,2001) 

WRITE (6,2002) (II, XXPRES( I 1.1) , I 1=1, NEQPOT) 
GO  TO  (1,2, 3, 4)  , ILOAO 

1 R£AC(5*2CG)XNXXA,0NXXA,XNXYA»ACCUR»RII 

120  'WRITE  (6, 511)  XNX  X A, ONXX A, XNXY  A, ACCUR 

OO  11  Il=l,N£QPOT 

11  AXPRES (II, 1 ) =X  X P RE  S ( II ,1) 

XNX=XNXXA 

DNX=DNXX A 

125  PF  AC  = 1 . 

' XN11=XNXXA 

GO  TO  5 

2 REAO(5,200)PFAC,OPFAC,XNXXA,XNXYA,ACCUR,RII 
WRITE (6,512)PFAC,0PFAC,XNXXA,XNXYA,ACCUR 

130  00  12  11  = 1, NEQPO  t 

12  AXPRES(I1,1)=PFAC*XXPRES(I1, 1) 

XNX=PFAC 

ONX  = J?F  AC 
X N 1 1=  PF  AC 

135  GO  TO  5 

3 REAJ( 5,200) XNXY A,DNXYA , X NX X A , A CCUR, RI I 
WRITE (6, 513) XNXYA.ONXY  A,XNXXA,A  CCUR 

00  13  11  = 1 ,Nf.QPOT 

13  AXPRES (II, 1)=XXPRES(I1 ,1) 

140  ‘ XNX=XNXYA 

DNX=DNXY  A 
PF AC= 1 • 

XN1 1=  XNXY  A 
GO  TO  5 

145  4 RE AD ( 5, 20QIPFAC, OP FAC, FNI»XNXY A, A CCURfRII 

- 00  14  11=1 jNEQPOT 

14  AXPRES (II, 1)=PF AC*XXPRES(I1, 1) 
XNXXA=FNI*AXPR£S (1,1) 

ONXXA=DPFAC*XNXXA 

150  XNX=PFAC 

DNX  = JPF  AC 
P1=AXPRES(1,1) 

0P1=0PFAC*P1 
X N 1 1-PF  AC 

155  WRITE (6, 51 4) PI, DPI, FNI.XNXX A, XNXY A, ACC UR 

5 IRk=R II 

IF(IRR.EQ.-Q)IRR  = 1 

XFNX=XNXXA 

XFN Y=XNXY  A 
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XFNP=?F AC 

TXNX=  13  0 0.*XNX 

RE  A 0 ( 5 , 1 GO ) NNN  » LNNN,  ILNW 

NWA VE=NNN 


***»,*»' 


♦ *»**V*CS****  * **♦*♦!  **♦*♦*«.*,♦♦***♦♦***♦*♦♦♦*♦*♦*****♦**♦**♦♦♦** 

' Wri  I T£  <b,505)NNN 

IF(LPRINT.NE.l) GO  TO  49 
W r<  1 T t ( b » 5 0 6 ) Cl » C 2 » C 3 , C4 * C5 » C 6 
WRITE(b,507)C7,C3,C9»ClQ»Cll»C12 
49  CONTINUE 
ILR=G 
LICON-1 
IPOTT  = 0 

CALL  SECONO(TIMl) 

WRITE  (6, 793) TIM1 
TIH2=TIM1 
TIM4=TIM2 
1 1 NN=  0 
555  LN=1 

CALL  AALOAD(ILOAO,NEQPOT,F  NI , XNX , XFNX, XFNY ,XFNP ,XXPRE S , AXPRES, M3, 
1114) 

IOET=IDOET 

CALL  POTERS (IO£T,NRHS, MAXN, AP, BP, CP, GP , PR , XP , CC  , M T , T1 , V 1, M AX  2 , 
1IXPM, UtTM, XFNX, XFNY,LN,NJ,NAW,N3W»NF) 

IF(LFRlNT.NE.l) GO  TO  101 
CALL  SECOND ( TI M3 > 

TIM1-TIM3-TIM2 
TIM2=TIM3 
TI  t-.4=  T I M 3 

WRITE (b, 201) NW A VE, XFNX  ,XFNY, XFNP.TIMl 

101  CALL  TRANSFtTl, MAXN, 1 ,LPRINT) 

IOET=IOOET 

CALL  AALOADJILO AO,NEQPOT ,F NI ,XNX, XFNX, XFNY ,XFNP .XXPRES, AXPRES,M3, 
1M4 ) 

I A MAX =1 
I8i1AX  = l 
A W N A X = C • 

3WHAX=0« 

I T£R=  G 

DO  102  J 1=  1 , K1 
A W M A X = A W M A X AW  M ( 1 , J 1 ) 

IFU1.EU.K1)  GO  TO  102 
8WMAX=3WMAX*3WM  (1,  Jl) 

102  CONTINUE 

DO  103  1 1 = 2, NEQPOT 
AW.1M=0. 

3WMM=0. 

30  104  J1=1,K1 
AWMM=AWMM«-AWMII 1 , Jl) 

IFlJl.tU.Kl)  GO  TO  1G4 
•3W  M 11=  d W M M §■  B W ‘1  < 1 1 , J 1 1 
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104  CONTINUE 

IF (ABSTAWMMI .LE.ABS(AWMAX) 1 GO  TO  32 
AW MAX  = AWMM 
I A :-1  A X = I 1 

32  IF  IABSMNMM)  .LE  . ARS(BWMAX)  I GO  TO  103 
UWMAX-BWMM 
I j A x = 1 1 
103  CONTINUE 
J A W M A X ■=  1 
J3  WMA  X= 1 

AWWH(1»=AWM(IAMAX,1I 
A A W W M = A W W M ( 1 ) 

BWWU<1>=3WM(IBMAX,1) 

0AWW.1-6WWM  <1) 

IF ( K1 • EQ. 1 ) GO  TO  1051 

DO  105  J1=2,K1 

AWWM( Jl»=AWM(IArtAX, Jll 

IF (A3STAWWM C J1 ) ) . LE . ABS ( A A WW M > ) GO  TO  31 
AAWWJ-i=AWWM  (Jl) 

JAWilAX- Jl 

31  IFCJ1 .EO.KIIGO  TO  105 
3WWM( J1>=BWN(IBMAX, J1J 

IF ( AdStBWWNt Jl)  ) .LE. ABS(BAWWM)  ) GO  TO  105 
iJAWW!4  = BWWMI  Jl) 

J BN  MA  X— J 1 

105  CONTINUE 
1U51  J C K A X = 1 

J O H A X = 1 

CCM  ( 1 l-CFM(IAMAX»l) 

DUU<1)=DFM(I3MAX,1) 

ACCN=CCM (1) 

A0D(1=  DON  ( 1 ) 

IF(K<2.EQ.l)GO  TO  333  ' 

DO  1 J 6 J 1=  2 » KK2 
CCM(Jl) =CFN( IAMAX»J1I 
DOK(Jl) =UFM(IBMAX,J1> 

IF (AbSICCM (Jl) ) .LE.ABS (ACCM) )GO  TO  30 
ACCM=CCM(J1> 

J C M A X = J 1 

30  IF(ABS(DDM (Jl) ) .LE.ABS (AODM) )GO  TO  106 
ADDN=  DDrl  ( Jl) 

JOMAX=Jl 
1C  6 CONTINUE 
333  LN=2 

ITER-ITERU 

IF  T ITER  .LE.10) GO  TO  113 

WRITc.  (6»  114)  ITER 

GO  TO  9999  - 

113  CALL  PCTERS(ID£T,NRHS,MAXN,AP,BP.CP»GP,PR,XP,CC,MT,T1, VI, MAX 2, 
II XP  H * D£  T M»  XFNX  * XFNY ,LN  » N J*  NAW , NBW ,NF> 

IF (LPRINT.NE.1) GO  TO  111 
CALL  SECOND (TI M3> 

TIM1=TIM3-TIM2 


' 


TIM2=TIH3 

WRITE  (6, 112)IT£R,NWA\/E,XFNX,XFNY,  XFNP,  TIM1 
111  CALL  TkANSF(T1,MAXN,1,LPR1NT> 

□0  115  J1=1,<1  ‘ 

IFIAWM(IAMAX.Jl) .NE.U)G0  TO  58 
AWCON(ITEP,Jl)  = 0. 

GO  TO  56 

5 8 AWCGN  (ITER,  J1J  = Ad5  UAWMII AKAX,  JD-AWWM  (Jl)  )/AWM  1 1 AM  AX,  Jl)  > 

56  IF { Jl .EQ.K1) GO  TO  115 

IF  16NM(I8MAX,J1)  .NE.  D.  )G0  TO  57 
3WCON(ITER»J1J-0« 

GO  TO  115 

57  CONTINUE 

dWCON  (ITER,  Jl)  - AOS  ( (3WM(  I UMAX,  Jl)  -3WWI-K  Jl)  )/BWM  ( I BM  AX , J 1)  I 

115  CONTINUE 
AWCH=AWCON(ITER, JAKMAXI 
3WCH=BWC0N(ITER, J3WMAXI 
I AWH= JAW MAX 
IOWH= JOWMAX 
00  116  J1=1,KK2 

IF  (CFM(IAMAX,J1)  .NE.O.TGO  TO  5*4 
CCONUTER,  Jl)=3. 

■GO  TO  53 

54  CCON  (ITER, J1)  = A0S( (CFM ( I AM AX , Jl ) - CCM ( J 1) ) /CFK( I AMAX, Jl )) 

53  IF  (DFM(I6MAX,J1>  .NE.O.  )G0  TO  52. 

000N( ITER, J1)=0  . 

GO  TO  116 

52  000N( ITER, J1)=ASS( (OFM (I 6MAX , Jl ) -DOM ( J 1 ) ) /OFM ( I SM A X , J 1 ) > 

116  CONTINUE 
CCH=CCON(ITER, JCMAXI 
ODH= J DON ( I TER,  J DMAX ) 

ICH  = JCI1AX 
IDH=JOMAX 

IF(LPRINT,NE,1) GO  TO  117 
WRITE (6, 1181  ITER, AWCH , BWCH ,CCH , JOH 
WRITE (6, 119) (Jl , AWCON ( ITER , J 1) , J1=1,<1> 

WRITE  (6, 119)  (Jl, 3WC9N  ( ITER,J1)  , Jl  = l,  KFOUR) 

WRITE (6,  119)  (J1,CC0N(ITEF, Jl) ,Jl=l,K<2> 

WRITE (6,119)  (J1,D00N(ITER,J1) ,Jl=l,KK2) 

117  IF  lAwCH.GT • ECON V ) GO  TO  194 
IF  (BWCH.GT.ECONV)GO  TO  194 
IFICCH.GT.ECONi/)  GO  TO  194 
IF(00H.GT.EC0NV)G0  TO  194 
IF (XNX.EC.XN11)  GO  TO  193 
IF  (QPNS70£T:1.GT  . J.)  GO  TO  193 

■ WRITE (6, 192) OETM.IXPM  , 

GO  TO  197  ' 

19  3 OPHS= DE  TM 
GO  TO  195 

194  IF l I T£R . LE . 2) GO  TO  196 

IF(  AWCONdTER,  IAWH)  , GT  . A WC  ON  ( I T£R-1 , 1 A WH ) ) GO  TO  197 
IF (3WC0N (ITER,  I BWH) ,GT .3WC0N  l ITER-1 , IBWH) ) GO  TO  197 
IF (CCGN(ITER,ICH) ,GT •CCONtITER-1, ICH) ) GO  TO  197 


SO 


IF  (JUUiMUlwR.IUri)  .GI  . uJJHll  fck-1,  lunn  00  10  1 92 

GO  FO  19b 

19 7 IF  (XNX.NE.XN11) GO  FO  198 
WRITE (6, 991) XNX 

GO  FO  9999 
19b  00  131  J 1= 1 » K1 

A W W N (Jl) = A W M ( I A H A X , J 1 ) 

IF  t Jl.EQ.KDGO  FO  131 
3WWM( Jl)  = BW,-|(I8MAX,J1) 

131  CONTINUE 

JO  132  J 1 = 1 , KK  2 

CCH (Jl) = CFM  ( 1AMAX.J1) 

132  0DN(J1)=0FM(IDMAX,J1) 

GO  FO  333 

195  CALL  AAL0AD< ILOAO,N£QPOT,FNI,XNX, XFNX, XFN Y , XFNP , XXPRES , AXPRES.M3, 
1M4) 

CALL  PO  F SN ( POT  * SIR  YU*  STRAU, STR YG, STR AG .XFNX.XFNY) 

CALL  SECONOC  TIM3) 

TIr!l=FIM3-TIM4 

tih2=tim3 

TIM4=FIM3 

NR  I Fc (6, 241) XFNX.XFNY, XFNP ,NWAV£, ITER, T INI 
WRITE  (6,242) PO T , STRYU, STRAU, STR YG, STR AG 
IF ( 10 ET  .EQ.l) WRITE (6,2  43 IDETM, I XPN 
’ IF (LMOU. N£ . 1)G0  FO  426 
CALL  TRANSF(T1,MAXN,2,3) 

426  CONTINUE 

IF(LNNN.EQ.2.ANO.LICON.N£.1G)GG  TO  556 
IFILICON.NL. 10) GO  TO  629 
ILR=1 
GO  TO  777 

b29  XNX 1=  XNX  ' 

IINN=I1NN*1 

IFCI1NN.LE.IRR) GO  TO  221 
WRITE (6,222) IINN 
GO  FO  9999 
221  CONTINUE 

XNX=XNX*DNX 

IF(TXNX.GT .XNX) GO  TO  244 
UNX=0NXX2. 

XNX=XNX-ONX 
ADN=ONX*iOO/XNX 
IINN= IINN-1 

IF(ADN.GT.ACCUR) GO  TO  244 
XNX=  X NX  1 
• GO  TO  819 
244  GO  TO  555 

198  IF(LICCN.NE.IO) GO  TO  429 
ILR=0 

GO  FO  777 

429  IF(LNNN.EQ.O)GO  TO  249 
IPOT  T = I POT  T* 1 
IF  (IPOFT .GT.l)  GO  TO  249 


S» 


r’U  i I - i-'O  T 
PXNX  = XNX-ONX 
249  AON=QNX* 100 ./XNX 

WRITt (6,545INWAVE,XNX 

TXNX=XNX 

XNX=X  NX  - DNX 

IINN-0 

IF  IAON.LE. ACCUR) GO  T 0 019 
DNX=CNX/2. 

XNX=XNX*ONX 
GO  TO  555 

819  CALL  AALOAO ( ILOAD, NEQPQT  ,f  NI , XN X,  XFN X , XFNY ,XFNP , XXPRES , AXPRES, M3 » 
1M4J 

WRITE 16, 785) NW A VE, XFNX, XFNY, XFNP, POT , S TRYU , STR AU, STRYG , STRAG 
C CALCULATION  OF  CRITICAL  WAVE  NUMBER 
566  WRITE (6* 20) 

IF (LNNN. EO.OJGO  TO  9999 
WRITE (6,584) 

ILR=1 

NW  PRI N=NWA  V£ 

UWAVc.  = NWAV£*-l 

IF(LNNN.EQ.2)POTT=POT 

P0TMIN=P0TT 

NMIN=NNN 

INWAVE=0 

ISTOP=0 

19=0 

POT=POTT 

777  CALL  AALOAOilLOAU, NEQPOT,FNI, XNX, XFMX, XFNY , XFNP , XXPRtS , AXPRES, 

1M3 , M4 ) 

IF (ILR.EQ. 1) GO  TO  778 

WRITEI6,582)NWAV£, XFNX, XFNY, XFNP  ' 

GO  TO  9999 

778  INWAVE=INWAVE*1 
1 1 NN=  0 

IF (LNNN. £Q  • 2 XNX  = XNX 

XNX=PXNX 

LICON=10 

IFIINWAVE.EQ.1)  GO  TO  391 
NWPRI  Im-NWAVE 
IFCI9.GE.DG0  TO  694 
IF  (POT MIN, LE. POT) GO  TO  139 
PO  TMlN=POT 
NM IU=  NW  A VE 
NWAVc=NWAVcH 
GO  TO  391 

139  IF(NWAVE.LE.NNN*l)GO  TO  549 
ISTOP=l 
GO'  TO  135 
549  19=19*1 

IF ( 19. Gf . 1) GO  TO  694 

NWAVc=NNN-l 

GO  TO  391 
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b94  IF (PCJMIN.LE.POT)GO  10  o95 
POTMIN=POT 
NM I N=  NW  A VE 
NWAVE-NWAVE-1 
GO  10  391 
695  1ST  0P  = 1 
GO  TO  185 

391  CALL  CGEFNNINWAVE) 

CALL  SECON 0 ( TI H2) 

Id  5 WRITE  (b,5311NWPRIN,PXNX,P0T, TIM2 
IF (NWAVE.LE. 0)GO  TO  9999 
IF  (ISTOP.NE.l)GO  TO  798 
WRITE  (6, 78  91  XNX , POTMIN  ,NMIN 
GO  TO  9999 

7 98  IFIINWAVE.GT .ILNWIGO  TO  9999 
GO  TO  555 

9999  R£AO(5»10QILCONTN 

IF1LC0NTN.E0.1IG0  TO  1111 
20  FORMATllHU 
10  FORMA  T ( 1H0  » 9A8  T 

60  FORMAT (//25H  BEGINNING  OF  NEXT  CASE  ///) 

100  FORMAT (1016) 

20  0 FORMAT  ( 6E12. 4) 

300  FORMAT  (//,2X, "NO.  OF  POI  NTS=*M  8 , 2X , "KFOUR=" , I 8 , 2X  , "BOUND.  CON  OF 
1 POINT  1 = “»I8»2X»" POUND. CON.  OF  POINT  NEQ POT  = “ , 18 >x 
LOO  F0RMAT(//,2X,"R=",E12.4,2Xf"XL=",El2.4,2X,"XH=",E12.4, 2X,"ELAS=", 
2E12. 4,2X,"XNI=",E12.4, 2X,"00=", £12.4, 2X,"EXXP=" ,E12.4) 

57  3 F0RMAT(//,2X,"XLAM0=”,£12.4, 2X , "YL AMU  = ", El  2. 4 , 2 X , "£X  = " , £12 . 4 ,2X , 
1"EY=",E12.4,2X, "RHGX=",E12.4,2X,"RH0Y=",E12.4) 

508  FORMAT!//, 2X»” THE  IMPERFECTION  SHAPE  IN  AXIAL  DIRECTION  IS"/) 

510  FORMAT  1//,2X, "THE  IMPERFECTION  FOR  CIRCUMFEREN TIAL  WAVE  “,I6/) 

52  0 FOR  M AT  (//,4X,”  POINT  ”,9X,  "LENGTH  ”,10X,"WZCOS",12X,"WZPCOS**,10X, 

1"WZPPCGS",13X,"WZSIN”, 1 0 X , "WZPS IN", 1 0 X , "W ZPPSI N ”/ ) 

509  FORMAT (110, 4E16. 6) 

5091  FORMAT (110, 7Elb. 6) 

500  FORMAT (//,2X,"DELTA=",E12.4,2X, "ALl=",tl2.4,2X,”GAl=",E12. 4,2X, 
2"AL2=",£12.4,2X,"BT2="»cl2.4,2X,"GA2=“,E12 .4) 

501  FGRMAT(//,2X,"H11=",E12.4,2X,"H12=",E12.4,2X»"H22=",£12.4,2X, 
1"Q11=",E12.4,2X  ,"Q12  = ",£1E.4 ,2X  ,"Q22=",E12 .4) 

502  FORMAr(//,2X,"011=",E12.4,2X,"O12=",E12.4,2X,”D22=",E12.4,2X,"OB2 
1=",E12.4,2X,"DB3=",E12.4,2X,"0&4=“,E12.4> 

50  3 FORMA  T(//,2X, “0L1=", hi  2.4, 2X,"OL2=", El  2.4, 2X,"DL3  = ", £12.4, 2X, 
1"0L4=",E12.4) 

50  4 FORMAT (//, 2X,"0A1=", £12. 4, 2X , "0 A2=", £1 2 . 4 , 2X ,"0 A3=" ,E1 2 . 4, 2X ,"0A4  = 
1" , £ 1 2 .4) 

112  FORMAT!//, 2X, "SOLUTION  FOR  I TER  AT  ION" , 18 , 2X, "N= ”,  18 , 2X , 
l"NXX=",E12.4,2X,"NXY=",El2.4,"FAC.LAr.L0A=",E12.4/2X, 

2"TIME  COMPUTATION  =", E 12 . 4 , 2 X , "SECON OS*V2 X ,"=  = = =================== 

3=r ================= ======= ===================== ====== ============= 

4=========================") 

118  FORMAT (//,2X,"ITER  = ", 1 8, 2X ,"AWCH=", E15 . 5, 2X ,"BWCH=", £15 .5 ♦ 
12X,"CCh=",£15.5 ,2X,"DDH=",£15.5/) 

119  F0RMAT(4 (I8,£15 .5)) 
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^>91  SOL  U 1 1 UN  IS  Nl)  I CONVtKGEU,  APPARtNI  L Y I HE 

1INIIIAL  LOAD  SHOULD  oE  LLSS  THAN",E12.4> 

2% 1 FORMAT l// , 2X»"N  ONL INE AR  SOLUTION  FOR  NXX=" ,E12 . 4 , 2X , "NX Y = " , E 12 . 4 , 
12 X, "FAC.LAT.cOA  = ",E 12. 4 , 2X  , " WA  V £=  ",  I 8 , 2X,"IS  OBTAINED  BY",  18  , 

22X ,”I TER AT  I ONS"/ 2X » "T I ME  COMPUTAT ION=",E12 .4,2X , "S ECON US"/ 2 X ," 

• sokokokckokokokokokokokokokokokukckokokokokokokokokokokokokokokokok 
40KOKUKOKOKOKOKOKOKO<OKOKOKC^OKOKOKOKQKOKOKOKOi<OKOKOKOK") 

192  FORMA  T ( // , 2X ,” T HE  DETERMINANT  CHANGES  THE  SIGN  DE T M = " , El 2 . 4 , 

12X , “I X?M  = " ,110) 

242  FORMAT (//,2X, "POTENTIAL  ENERGY  = ",  El  7 . 7 ,2X . "END  SHORTENING  U FOR 
1Y=0.  IS", E15. 5, 2X, "AVERAGE  END  SriOR.  U = ",  E 15 . 5 / 2X  , "END  SHORT, 

2 GAMA  FOR  Y=0.  I S" , El  5 . 5 , 2 X, "A V ER AGE  END  SHORT.  G AM A = " , El  5 . 5 ) 

722  FORMAT (//,2X,"END  OF  THIS  CASE  BECAUSE  NUMBER  OF  LOAD  POINTS 
1 IS  GREATER  THAN", 181 

545  FORMAT (///»2X»"FOR  WAVE  NUMBER" , 1 8, 2X, "THE  LOAD",  £12. 4,2X ,"IS 
1 OVER  THE  LIMIT  POINT") 

785  FORMAT (//,2X, "CRITICAL  LOAO  FOR  WAVE  = " , 18 , 2X , " ARE  N XX= " , El  2 . 4, 2X , 
1"NXY=",E12.4,2X ,"L AT. LOAO=", El 2. 4/2X, "POTENTIAL  ENERGY="»E17.7»2X, 
2"END.SHOR.  U FOR  Y = 0.  I S", El  5 . 5 , 2 X , “ A VERAG E END.  SHOR.  U=",El5.5/ 
32X  ,"cNC. SHOR.  GAMA  FOR  Y = 0.  I S"  , E 15 . 5, 2X , " AVER A GE . END . SHOR . G AMA=" , 
4Z15.5/2X ,"OKUKOKOKOKOKOKOKOKOKOKOKO<UKOKOKOKOKOKOKOKOKOKOKOKOKOKO< 
50KCK0K0K") 

58  h FORMAT  t//,2X, "DETERMINATION  OF  THE  CRITICAL  CE RCU MFERE NT  I AL  WAVE 
1 NUMBER  ( N)  "/// 2 X,  ••***♦*****♦**♦***♦*  + ♦**♦**♦♦♦***♦***  ,******♦**** 

58  2 FORMAT {//, 2X,"FOR  W AVE=" , 1 8 , 2X , "NXX  = " , E12 . 4, 2X , "NX Y=", E12. 4, 2X, 
1"LAT. L0AU=",E12 . 4 ,2X,"TH£  SOLUTION  IS  NOT  CONV  E RG£0"/2  X ," 

2PRU  3 A 3L  Y THE  CRITICAL  LOAD  IS  SMALLER  THAN  THESE  LOADS") 

581  FORMA T (//,2X, "WAVE  NO . =" , 1 10 , 2 X , "LOA D= ", E 1 5. 5 , 2 X , "P OTE N TI AL  ENERGY 
1=" , El  7. 7 , 2X , "EL A3SE D T IME  = ", El 2 . 4 , 2X , "SECONDS") 

789  FORMAT (//,2X,"F0R  LOAD", E12. 4, 2X,"THE  MINIMUM  POTENTIAL  ENERGY 

1 IS", E17.7,2X,"AND  THE  CRITICAL  WAVE  NUMBER  IS  ",I6) 

511  FORMAT <//, 2X, "THE  INITIAL  AXIAL  COMPRESSION  IS  NX X = ", £ 1 2 . 4 , 2 X , 

1”T HE  INCREMENT  IS  DNXX  = " , E 1 2 . 4/ 2X , "T HE  LATERAL  LOAD  IS  FIXED  AND 

2 THE  PERMANENT  TORSION  IS",E12 . 4/2X, "THE  ACCURACY  = ", El 2 . 4, 2X  , 
3“PcRCENTS'V) 

512  FORMAT!//, 2X, "THE  FACTOR  OF  INITIAL  LATERAL  LOAD  IS" , E 12 . 4 , 2 X ," 
1THE  FACTOR  OF  INCREMENT  I S", £ 1 2 . 4/ 2X , "THE  FIXED  AXIAL  COMPRESSION 
21S", £ 12 . 4 , "A ND  THE  FIXED  TORSION  IS" ,E 12. 4/2X , "THE  ACCURACY  = ", 
3E12. 4, 2X, "PERCENT"/) 

513  FORMAT J//.2X, "THE  INITIAL  TORSION  =",E12. 4/2X, "THE  TORSION  INCREME 
1NT  =",E12.4/2X, “THE  FIXED  AXIAL  L OAD  = ", El  2 . 4/2 X , "THE  ACCURACY 
2E12. 4, 2X, "PERCENTS"/) 

5.14  FORMAT  !//,2X, "THE  LATERAL  LOAO  IS",E12. 4, 2X,  "THE  INCREMENT  IS", E 
112. 4/ 2X, "THE  AXIAL  COMPRESSION  AND  THE  LATERAL  LOAD  ARE 

2 CONNECTED  BY  THE  F AC T OR”, El 2 . 4/2 X , "THE  INITIAL  AXIAL  COMPRESSION 

3 IS",E12.4,2X,"THE  FIXED  TORSION  IS" , E 12. 4/2X , " THE  ACCURACY  IS", 
4E12. 4, 2X, "PERCENT"/) 

505  FORMAT (//,2X, "THE  CIRCUMFERENTIAL  WAVE  NUMBER=",I6) 

506  FORMAT ( // , 2 X , "C 1= ", E 1 2 .4,2 X, "C2 =" , El 2 . 4, 2X ,"C3=",E12.4 ,2X,"C4=", 
1£12.4,2X,"C5=“, E12.4,2X,"C5  6=",E12.4) 

507  FORMA T(//,2X,"C7=",E 12.4* 2X, "C8=" »E12 . 4, 2X ,"C9=", El2. 4 . 2X , 
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l"ClU-"»t-12.4,2X,”Cll  = ”»E12.4,2X»C12=  *£12»4JI 
793  FORMAT  ( // » 2X »"E L APSEU  TI  !1E=”  , £ 1 2. 4, 2X>  "StCONOi"/) 

11  €*  FORMAT  <//,2X,,4£NO  OF  THIS  CASE  BECAUSE  ITER  CRt  A TER  THAN  ”i 

24  3 FORMA  T ( // , 2X ,”DE  TERM  IN AN  T= “ » £1 2 • 4 » 2X ,"IXPM=" , 1 1 0) 

2C1  FORMAT  <//,2X,  "INITIAL  SOLUTION  (FROM  LINEAR)  FUR  N=”,I8,2X, 
1"NXX=",E12 .4 ,2X ,"NXy  = ",E12.4 ,2X,“LAT . L OAO  = ",  El  2 . 4/ 2 X , " T IME 
2 C JMP U T A riON="»  E 12 . 4 » 2X,  "SEC0NQ5*V2X , ”= ==  = -==  = = ======  ==■==  = = 


2Dtll^F0kMATT//*2X»"  POINT  "»9X,"P  RE SSURE"/) 
2502  FORMAT (110 , E12,  41 
END 
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SUBROUTINE  A ALC  AtM  iLC-TOt  NtQPOI  • FNI  t XNX  » XFNXf  XFNYt  XFNPf 
. lXXPRtS,  AXPRES.'il  . M 

JIMENSION  XXPKESCilt  'O  . AXPKtS ( Ml ,M2» 

GO  TO  (1,2, 3,4)  , IlOAO 
1.  XFNX=XNX 
GO  TO  5 

2 XFN°- XNX 

UO  6 11  = 1 ,NEQP0T  _.,,T 

6 AX  PRES  (II,  1)  =XF  \’P*XXP.\Li  ( 1 1 » 1) 

GO  TO  5 

3 X F N Y = X N X 
GO  TO  5 

4 XFNP=XNX 

DO  7 11  = 1, NEOPO  T „.,#r4  .. 

7 AXFRES(Il,l»=XFNP*XXPRtS (11,1) 


XFNX=FNI*AXPR£S(1,1) 
5 CONTINUE 
RETURN 
ENO 
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FUNCTION  AE3(I, J.L.U. JP.N2.N3.N4,  LL> 

C • N 4- 1 FOR  0(JP,I+J)  Oft  BIJP.I-JJ 
C N4=2  FOR  d(JP,J) 

C LL=1  FOR  A LL=2  FOR  B.C.D 

CONMON/FOUR1/KFOUR.K1, K2.K3, K4, K5 , K6 , K7 , K 8 » K9 
CUMNCN/F0UR2/K1 0.K11. K12 , K13 , K1 4 , R15 , KK2 
DIMENSION  B(N2,N3) 

IF (L. GT. 3)  GO  TO  10 
I1=I*J 
12  = 11 
uO  TO  23 

10  Il=IAdS ( I-J) 

12  = 11 

20  IF (N4 .EQ.ll  60  TO  120 
12  = J 

GO  TO  100 

12 3 IF(Il.LE.KFOUR)  GO  TO  100 
AL 3=0  . 

RETURN 

100  LF (L . LE . 3)  GO  TO  110 
ETA=1. 

IF(I.cQ.J)  cTA=0. 

110  GO  TO  (11, 12, 13, 14. 15, 16. 22. 23, 24), L 

11  R1=I1**2 
GO  TO  17 

12  R1=J**2 
GO  TO  17 

13  R1=2.*I1*J 
GO  TO  17 

14  Rl= ( 2 . -E  TA) *11**2 
GO  TO  17 

15  ftl=(2 .-ETA) *J**2 
GO  TO  17 

16  IF<I-J.LT.G>£TA=-1. 

Rl=-2.*ETA*J*Ii 

GO  TO  17 

22  IF(I-J.LT.G)£TA=-1. 

Rl=-cTA*Il**2 

GO  TO  17 

23  IF(I-J.LT.0)£TA  = -1.~  

R1=-£TA*J**2 

GO  TO  17 

24  R1=2.M2.-ETA)*I1*J 

17  IF(LL. £0.1)12=1241 
AL8=Kl4t3  (JP,  12) 

RETURN 

END  i 
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SUBROUTINE  C«J£FNN  (NNN) 

COMMuN/GEOM/RR, 00,H11, H12,H22,Qll,Cl2,Q22, D11.D12.U22 
COi.;-13N/FACTOrJ/Cl,C2,C3,C4,C5,(;e>»G7.Ctt,C9.ClO»ClltCl2 
Cl  = »iNN/RR 
C2=r;i**2 
C3=G2**2 
C4  = Cc  *Ci  1 1/011 
C5  = C2/UR*011I 
C6 =2. *0J*H12*C2 
C7=2.*Q12*C2 
Cfl=JD*H22*C3 
C9-C3 / ( *♦  • *01 II 
CIO  =Ut2  * C3 
C11=2.*012*C2 
C12=022*C3 
RETURN 
ENO 
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SUBROUTINE  COEFF {£X,£Y ,XLAMO,YL  AMD,RHQX,RHOY,ELAS) 

Coni', ON/FOURt/KF OUR,  Kl»  K2,X3,K4,  <5,i<6,K7,<3  ,K9 
COMMQN/FOUR2/K10,<11«K12,K13»K14»K15»KK2 
COMNON/GEOM/RR , OD, Ml 1 , HI 2 , H22 , J 11 , Ql 2, Q22 , 01 1 , 0 12 , D22 
CQMM0N/FACT2/0L1 ,0L2,0L3 ,CL4 , 0 A 1 , 0 A2 ,0  A3 , 0 A4 , 0 J 2 , Ob3 , OB  4, XNI , EX XP 
COMMON/C  INI G/Nt  uPOT  *NI (5  0 & ) 

COMilON/F  IOFR/OE  L T A , AL1  ,GAl,AL2,BT2,GA2 

COMMON/F ACT  3/XL , XH 

KK2  = 2 .*KFOUR 

K1  = KF  OUR  + 1 

K2=Ki +1 

K3=2.*KF0UR+1 

K4=K3+1 

<5=4*KF0UR+1 

K6=K5+1 

K7  = 6*  KFOUR+1 

<3  = K7  + 1 

<9=  7* KFOUR+2 

Klu=K9+l 

Kll=J»KF0UR+2 

<12=<11+1 

Kl3=lU*<F0UR+2 

K14=<13+1 

K15=12*<F0UR+2 

XN2=XNI*#2 

XH2=XH**2 

0ALFA=(1.+XLAM0) * ( 1. +YLAMO) -XN2 
DQ=ELAS*XH+*3/ ( 1 2. * ( 1 . -XN2 ) > 

EXXP=ELAS*XH/tl .-XN2) 

HI  1=1  . +RHOX+  12.  *EX**2*XLAMD*  11. +YLAM D-XN2 ) / (XH2+DALFA) 

H 2 2=1. +R HOY +12. *EY ++2+YLAMQ* (1. +XLAMD-XN2) /(XH 2*0 ALFA) 

HI  2 = 1 .+12. *XXI*EX*£Y+XLAMD*YLAMD/ (XH2+DALFA) 

G11=-XNI*EX*XLAH0/0ALFA 

(12  2 = *XNI*£Y*YL  A HO/DALFA 

Ql 2= ( (1.+YLAND1  *EX*XLAMD+  { 1.  +XLAMO) *cY*YL  AMO) / C 2.+UALFA) 

011  = (1  . + XLAHD) / (OALFA+EXXP) 

012=  < (l.+XLAHO)  * (l.+YLAMU)-XNI)  / < OALFA* <1 . -XNI ) *EXXP) 

022= ( 1. +YLAMO)/ (OALFA+EXXP) 

0L1=0D*H11 

0L2=00*XNI* (1.*  (EX*EY*XLAhU*YLAM0*12.) / (XH2+DALFA)  ) 

0l3=-£X*XLAHD* ( l.+YLAMO) /OALFA 
DLh=XNI*£X*XLAM D/0 ALFA 
DA  1=  ( l.+YLAMO) /(DALFA*£XXP) 

0A2=-XNI/CDALFA*lXXP) 

OA  3=  (l.  + YLA:iQ)*EX*XLAMD/DALFA 
DA4  = -XNI*EY*YLAI'I0/UALFA 
D32=( l.+XLAHO) ✓ (OALFA+EXXP) 

OB3=-XNI*EX*XLAHO/OALFA 

094= 11. + XL AMD) *£ Y+YLAMO/UALFA 

ni(l)=2*K15 

MI ( Nc  QPOT) = 2*K 15 

N£Q1=N£QP07-1 

DC  10  11=2, NEQ1 


Ml  ( ID— K15 
10  .CONTINUE 

OE ».  TA  = XL/  ( NE  QPO  F — 1 > 
AL1  = -1./12.*JELTA) 
GA 1 = 1 ./ < 2.*0ELT  A) 

AL 2= 1 . / ( QELT  A**  2 > 

3T2  = -2./OEl.TA**2 

GA2=l./OELTA**2 

RETURN 

ENO 
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SUOROUT INE  dOUNDR  (dS,  OT  , faG , IN  , XNXX,  XNXY,  L S,  Ml , NJ  , NA  W,  NOW , NF  * 

. 1LN.NRHS) 

G0MMGN/F0UR1/KF0UR,K1,K2,K3,K4, K5,K6,K7,KB  ,K9 

CUMMON/FOUk2/K10,Kll,K12,K13,Kl4,K15,KK2 

COM. ICN/G  EOM/KR,  DJ,  HI  1 , HI  2 » H22  , Ul  1 . 01 2 , Q22 , 0.1 1 , 0 12 , 022 

COMMON/FACTOk/C 1,C2,C3,C4.G5.C6  , C 7,  C8 , C9, C 1 0, C 1 1,C12 

COMMON/F  ACT2/DL 1 ,DL2 , DL3  » 0L4 , DAI , DA2 , 0 A3, 0 A4, Jd2 , 0H3 , 004, X NI ,£XXP 

C JMMQ N/PRE  Sl/A W M (ICO*  31 , AW HP (10  0. 3) , AWKPP( 10  0,3) 

COMMON/ PRES2/BWM ( 100, 2) , OWMP (10  0, 2) , 8WMPP ( 10  0, 2) 

COMMON/PRE S3/CFN (100,4), CFMP (10  3, 4) , CFMP P ( 100, 4) 

COMHON/PRES4/DFM (100,4) ,OFM? (100, 4)  , OF MPP  (100,4) 

COMMON/PP.ES5/AWZ  (10  0, 3)  , AWZP l lu  0 , 3) , AwZPP  ( 100, 3) 

C OMMG N/ PRE S6/B W Z (100,2),  BWZP ( 10  0, 2)  , BW  ZPP  ( 10  0 , 2 ) 

DIMENSION  8S (Ml , Ml ), 8T  (Ml , Ml ) , 8 G ( Ml , NRHSI 
Al=l. 

IF ( LN . EQ. 1 ) AL=  0 « 

IF (LS.EQ.ll) LS=  3 
DO  1 11=1, K15 
3G(  II  ,NRHS)  = 0. 

DO  1 J1=1,K15 
3S(I1,J1)=G. 

3T ( II , J1)=0. 

1 CONTINUE 

IF(LS.NE.IO)  GO  TO  100  * 

DO  2 1 1=  1 , K3 

2 BS(I1,I1)=1. 

BS<<t5  ,K3)=DL 4*Q11/D11-0L1 
3G (<6,NRHS)=XNXX+AWZP< IN,1) 

E4  = -C2/2.*Di_  4/011 
E5=-C2/2  • 

DO  3 J=1 , KFOUR 

JS=J**2 

J1=J+1 

BT(K8 ,J1)=E4*JS* AWZP (IN, Jl) 

3T (Ko,K1+J)  = £4*JS*eJWZP(IN*J) 
or  (kd,K3+J)=E5+  J 3* AWZP  (IN,  Jl) 

3T (K3 ,K5 +J)=E5+JS*JWZP(IN, J) 

3 CONTINUE 
K81=K8+1 

1=0  

DO  4 I1=K81,K9  • 

1=1+1 

OS ( I 1 ,K8+I)=-DL1 
3S( II ,K11+I)=-0L4 

t3T  (I1,K1+I)  = XNXY+C1  + I ' 

JG ( II ,NRHS)=XNXX*AWZP( IN.I+l) -XNXY*BWZ( IN, I) +C1+I 
30  5 J=1,KK2 

ar ( Il,K3+J)=3T( I1.K3  + J) +E5+ ( ALE (I ,J,2, AWZP , IN,N J,NAW,1 , 1)  ♦ 

1AL  0(1  ,J,  5,  AWZP,  I N,  N J,  N AW  , 1 , 1 ) ) 

OT (II  , K5+J)=OT (I1,K5+J)+25*(AL3(I,J,2,BWZP,IN,NJ,N3W,1,2)+ 
lALd (I , J ,3, BWZP, IN,HJ,NOW,l,2)> 

5 CONTINUE 

4 CONTINUE 


* 
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I=* 

DO  6 1 1 = <1 0 , K1 1 

• 1*1*1 

i3S(Il,K9*I)  = -DE  1 
GS<I1,K13H>=-QL4 
• BT  dl , 1*1) =-XNXY*Cl*I 

BG  < 11  ,N?H5> =XNXX*BHZP( IN, I)  +XNXY*C1* I* AWZ ( IN, 1+ 1» 

00  7 J- 1 » KK2 

dT (II ,K3+J)=oT ( I1,K3  + J)*£5M  AL3  (I , J,2,  BWZP,IN, NJ,N3M,1,2>- 
1ALB (I ,J, 5, EM ZP, IN,NJ,N3H,1,2) J 
UT  (I1,K5*J>=GT(I1,K5*J)*£5M-ALB(I,J,2,AWZP,IN,NJ,NAM,1,1)* 
lALBd,J,5, AMZP,  IN,NJ,NAM,1,1>> 

7 CONTINUE 
6 CONTINUE 

00  ft  I1=K4,K7 

8 BS ( 1 1 , 1 1 i = 1 • 

00  9 1=1, KK2 
BS(K11*I  ,Klld>  =DB2 
3S(K1 3*1 ,K13+II  =082 

IF  (I.  siT.  KFOUR)  GO  TO  55 
3S  ( K 1 1 * I » K ti*  I)  = 0B3 
3S(K1B*I,K9»-I)  = D33 
55  00  10  J=1,KF0UR 

3T(K>ll*I,J*l)  = BT(Kll*I,J*l)*£5*(ALB(I,J,2,AWZP,IN,NJfNAM,ltl)* 
1ALJCI .J.5, AMZP,  IN,NJ,NAM,1,1) ) 

BT  CK11+I ,Kl*J)  = DT(Kll+ItKl*J) *E5* < ALB < I , J , 2, DWZP, IN,NJ,NdM,l ,2) ♦ 
lAl_d(I  ,J,8,BWZP,  I iT,  NJ  , N BH  , 1 , 2 ) ) 

dT( K1 3*1, J*1)=BT(K 13*1, J*li*E5* ( AL6( I , J ,2 , BWZP , IN , N J,  N BM , 1 , 2 ) - 
1AL6(I,J,8,BMZP, IN,NJ,NBM,1 ,21) 

3T(<12*I  ,K1*J) =6T(K13*I,K1*J>*£5* l -ALB (I, J ,2 , AM ZP » IN, N J , NA W, 1,1)* 
1AL.KI  , J,5,  AMZP,  I N,  N J , NAN  ,1,1)) 

10  CONTINUE 
9 CONTINUE 
RETURN 

100  IFILS.NE.9)  GO  TO  200 
E<,=0Ll-JL4*m/011 
E5=UL4/ (RR*011) 

£6  = 0L4*C2/ (‘♦.•Dili 

BT ( l , Kd ) =£4 

BS (K8 ,K3)=E4 

BT i 1 , 1) = -Eb 

dS (Kd , 1) =-E5*XNXX 

dG (Kd ,NPHS)=-XNXX*AWZP (IN, 1) 

DO  11  J = l,K.FOUR 
• JS  = J * *2 
J1=J*1 

BT  (1,  J1»=BT  11,  J1)*2.*E6*JSMAL*  AWM(IN,  Jl)  ♦ AWZ  (IN,  Jl)  > 

BT (1, K1*J)=BT(1 ,K1*J) *2.*tb*JS* (AL*dNM (IN, J)*3WZ( IN, J) > 
dG  ( 1»NRHS)=i3u(  1 , NRHS)  *E6*  JS*  AL*  ( A WN  ( I N , Jl ) **2+0 MM  (IN,J)**2) 
3S(K8,Jl)=bS (K  8 » J 1 ) t2»*£6*JS*(AL*AWM(IN,Jl)+AHZ ( I N , J 1 ) ) 

BT  ( Kd  , J 1 ) = BT  (Kd  , J 1 1 ♦ 2 . * Eb*^SM  AL*  AWN P dN,  Jl)  *ANZP ( IN, Jl) > 
0SUd,Kl*J)=GS<K8,Kl*J)*2.*£6*JS*  ( AL  *BWM  ( I N,  J)  * BM Z < IN,  J > > 

BT  ( K& , K 1 * J ) =BT  (Kd,Kl+J)*2.*E6*JS*(AL*8MMP(IN,J)  *BNZP(IN,J)  I 
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CG  U3,NRHS)=.UG  l K8,NRHS)f2.*E6+JS*AL*  < A WMP  < IN,  J 1 > * A WM  U N , J 1 ) ♦ 
• liiWMPl  IN,  J)  *DWM<  IN,J>) 

11  CONTINUE 

00  12  I = 1 , KF OUR 
11=1+1 

31  (Il,K8+I)=0Ll 

nr ( II , Il)=-DL2+C2*I+*2 

OT < II ,K11+II=DL4 

JT  IK1  + 1 ,K9+I)=DL1 

bT (<1+I,K1+I)=-DL2*C2*I**2 

or  (Kl  + I , K13+I)  = UL4 

US lXd*I,K8+I>=0Ll 

□ S(<o»I,Il)  = -C2  + I**2MDL2+2.*DDMl.-XNI))+XNXX 
US(i<8+I,Kll  + I)  = Dl.4 
3TU6  + I,K1  + I>  = -XNXY  + I*C1 

3G(K8+I,NRWS)=-XNXX*AWZP  (IN,  11)  + X NXY  *1  *C1  * BWZ(  I N,  I) 
BS(K9+I,K9+I)=DL1 

0S(X9  + I,K1  + I)  = -C2*I**2M0L2  + 2.*DD*(1.-XNI>  J+XNXX 
3S(<9  + I»K13+I)  = 0L4 
BT U9  + I,I1)=XNXY*C1*I 

uG (K9+I,NRHS)=-XNXX*BWZP(IN,I) - XNX Y*C 1 * I* A WZ ( IN  , 1 1 1 

12  CONTINUE 

□0  13  1=1, KK2 
BI  (<3+I,K3+I)=l. 

OS (Kll+I ,K3+I)=1. 

3T  <K5+I,K5+I>=1. 

BS1K13+I  ,K5+I>  = 1. 

13  CONTINUE 
RETURN 

200  DO  14  11=1, <3 

14  dT ( II  , Il)  = l. 

IFU.S.CE.4)  GO  TO  300 
J=  0 . 

DO  15  I1=K8,K11 
J=J  + 1 

15  35(11,0)=!. 

GC  TO  400 

300  BT ( K3 , Kd) =1 . 

DO  16  I=l,*FOUR 

ii=i+i  ' : 

3T  U3  + I, K8+I)=OLl 
3T  ( K 3 + 1 » KU  + I)  = 0L4 

3T (K J+I, K9+I)=DL1  ’ * . 

BT  (K9  + I,  K13  + I)  = 0*_4 
1F(LS.EQ.1.0R.LS.EQ.3)  GO  TO  16 
JT (K3+I»K3+I)=-DL3*I**2*C2 
BT  ( K9 + I • K5+I ) =- OL3* I ** 2*C2 

16  CONTINUt 
400  E4=-C2/2. 

00  40  1=1, KK2  ' 

IS=I+»2 

11=1+1 

GO  TO (21,22,23, 24,25,26,27,28) ,LS 
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21  US ( < 3 ♦ I , <3  + 1 1 = 1 « 

. US(K3+I»K5+II=1, 

UTUll+I, <3+11  = 1. 

3TU13+I,K5  + I)  = 1. 

GO  TO  4 0 

22  US (<J  +1 , K3  + 1 1 = 1 . 

8SU5  + I ,<5+11=1  . 
nS(Kll+I  ,<11+11  =082 
US  (K1  3 + I ,<13+11  =082 

IF ( I • GT • KF OUO  GO  TO  41 
dS(Kll+I  ,K8  + I)  = DU3 
US  (K13+I  ,K9+I1  = DB3 

BS(K11+I,I+1)=-Q34*IS*C2+1./RR  ' 

3S  ( <13*1  ,K1+I)=-034*IS+C2+1./RR 

3S (Kll+I, 11=2. *E4*IS*AMZ(IN, 111 + 6S (< 1 1 + I , 1 I 

3S(K13+I,1)=BS  (K13+I,1)+2.*E4*IS*CMZ(IN,  I) 

41  JO  42  J=  1 » KFOUR 
J1=J+1 

US  (<11  + I , J1)=BS  (<114-1,  J11+E4+  ( ALB  (I»  J,  1,  AWZ,IN  , NJ,  NAM,  1,11  ♦ 
1ALB(I  , J,  4,  AWZ,  IN,UJ,NAW,1,  1)  1 

US  (<11  + I ,K1  + J>  =35 (Kll+I,  Kl+J)  + E4*  (ALB (I  ,J  , 1,BMZ  ,IN,NJ,N3W,  1,21* 
1AL3(I  ,J,  7,  BWZ,  IN,NJ,NBW,  1,2)  1 

oS (Klo+I,Kl+J»=BS(K13+I,Kl+J)+E4* (- ALB ( I,  J , 1 , A WZ, IN , NJ , NAM , 1 , 1 1 ♦ 
1ALU(I,J,4,AMZ,IN,NJ,NAW,1,11) 
8S(K13+I»Jl)=8S(K13+I»Jl)*£4+(AL3(I»Jfl»BWZ»IN»NJ»NBMfl»2)- 
1AL3 (I , J , 7 , 3MZ, IN, NJ.NBW, 1,2)1 
CONTINUE 
GO  TO  4 j 

3T  (K3+I, Kll  + I)  = D62 
UT  (K5  + I, <13+11  = 052 
IFII.GT. KFOUR)  GO  TO  43 
3T (K3+I»K8+Il=0B3 
BT  {<5  + 1 »K9  + I)  = 053 
6T  (Kll  + I ,K3+I)  = 1. 

BT  {<13  + 1 » <5+  I)  = 1 . 

GO  TO  40 

BT(K3  + I, <11  + 11  = 082 
BT (K5+I,  <13+11  = 032 
3T (<J+I,K3+I)=-0A2*IS*C2 
3T  (<p+I,<5+I)  = -0A2*IS*C2 

bS(<ll+I ,KJ  + I)  = -IS*C2* (DA2  + 2./( (l.-XNI) +EXXP)  1 
8S (Kll+I ,K11+I) =332 

9S (K13+I ,K5+I)=-IS*C2* (OA2+2 .✓ ( (1 .-XNI 1 *EXXP)  1 
3S ( < 1 3 ♦ I , <13+11 =0B2 
IFII.GT. KFOUR)  GO  TO  44 
3T (<5+I, <8+11=033 
3T (<5+I»K9+I)=UB3 
BS  (Kll  + I ,<8+11  = 033 
BS (Kll+I , I1)=-0B4*IS*C2+1./RR 
dS(<13+I, <3+11=083 
BS (Klo  + I ,K1+ 11= -DB4+IS+C2+1./RR 
US(<11+I, 11= dS (Kll+I, 11+2. *E4*IS*AWZ(IN, 111 
aS(K13+I,l)=BS(K13+I,l)+2.*E4*IS*3WZ(lN,Il 


42 

23 

I * 

43 

24 


44  30  45  J= 1 « KF OUR 
J1  = J + 1 

JS  <K11  + I,J1)=3S  (<11+I»J1)+E4*  ( A L b ( I * J < 1»AWZ»IiM»  N J » N A W t 1 1 1 ) * 

1AL13(  I t J,4,  AWZ.  IN,NJ,NAW,  1,  D) 

JSUll  + I»Kl  + J)=3S(Kll  + I,Kl+J>  + E4*(ALd(I,J,l,BWZ,IN,NJ,N3W,l,2)  + 
1AL3(I , J,7,dHZ,IN,NJ,NBW, 1,2) I 

dS(Kl3  + I,Kl  + J)  = 8S<i<13+I.Kl+J)+E4*  (- Ai_B  1 1,  J , 1 , A WZ,  I N t N J , NA  W t 1 1 1)  + 
1ALB(I,J,4,  AWZ.IN.NJ,  NAW.l,  D) 

3S  <K1  3*1 » Jl)  =3S  ( K13+I,  Jl)  +E4 * < ALB ( I , J , 1 ,BW Z, IN , NJ , NBW, 1 ,2) - 
1ALi3  II , J,7,  bWZ»  IN,NJ,  NBW,  1,2)  I 

45  CONTINUE 
GO  TO  40 

25  9S(K3+I,K3+I)=1. 

3SU5+ItK5  + IJ  = l.  . 

BT(K11+I,K3+I>=1. 

8T (K13+I ,K5+I)  =1. 

GO  TO  40 

26  BS(<d+I ,<3+I)=l. 

8S{<5+I,K5+II=1. 

BS  (Kll+I  »K11H)  =082 
JS (<13  + I .K13  + II  =032 
IF  (I.  GT  .KFO'JR)  GO  TO  40 
’3S(K11  + I»K8+I)  = 333 
dS(K13  + I,K9+I)  = 083 
GO  TO  40 

22  dT<K3H,K2  + I>  = l. 

BT(K5+I,K5+I>=1. 

3T  (<11  + I ,K11+I>  =082 
dT  (<1  3+1  t K1 3 + 1 > =032 
IF (I. GT.KFOUR)  GO  TO  40 
3TIKll  + I»K8  + II=l)B3 
3T  (K13  + I ,K9  + I)  =033 
GO  TO  40 

28  BT <K3+I,<11+I)=DB2 
9T  (K5+I,  K13  + I>  = 032 
3T  <K3  + I»  X3  + I)  = -3A2*IS*C2 
dT (K5+I , K5+I)=-0A2*IS*C2 

3S<KllH,K3+I)=-IS+C2*<0A2  + 2./((l.-XNI)*EXXP>) 
dS(Ki.l  + I»Kll+I)  =032 

OS<<lJ  + I,K5+I)=-lS  + C2MDA2  + 2./((l.-XNI>+£XXP)) 

BS  (K1 3 + 1 ,K13+I) =332 
IF  (I . GT.KFOUR)  GO  TO  40 
dT ( < 3 + 1 , K8+I)=033 
3T  (K5  + I , K9+I)=D33 
OS  ( Kl  1 + I i K 8 + I)  = 0 33 
3S(Kl 3+1 » K9+  I) =Od3 
40  CONTINUE 
RETURN 
£NO 


SUBROUTINE  RSTGtRtStT  , G , J P , XNX X, XNX Y , Ml , N J, NA W , NOW , NF , LN . NRHS) 
• GOKNON/F  OUK1/KF OUR , < 1 , K2 * K3 , K4 , K5 , K6 * K 7 , K 8 ,K9 
COMMON/FOUKZ/KIO ,K11  ,K12 ,K 1 3 ,K1 4, K15 , KK2 
CO.-MCN/PRES1/AKM  (100,3)  , AWMP  (10  0 , 3) , AWMPP  (10  0,3) 

• OOHMGN/PRSS2/L.WH  ( 1 0 0 , 2 ) , BWKP  ( 1 u 0 , 2) , 3WMPP  (10  0,2) 
CO'1MUN/PRtS3/CFM(lC0,4),CFMP(100, 4>,CFMPP  (10  0,4) 

COMHO N/PRCS4/0F H ( 1Q0 » 4) » UFMP ( lu  0 , 4) , uFMPP (100,4) 
C0HM0N/PRES5/AWZ  <100, 3) , AW ZP ( 1 0 0 , 3) , AW ZPP ( 10 0 , 3 ) 
GUMM0N/PRES6/3WZ ( 1 0 0 , 2 ) , 3UZP ( 1 0 0 , 2) , dWZPP ( 100 , 2 ) 
30l'li10N/XX(.OAU/AXf>RES(100,3),&XPRtS(100,2) 

COMMON/ GEOM/RR,  OO, HU, Hl2,H22, 011,012,Q2  2,011, 012,022 
CUMH3N/FACT0R/C1 , C2 , C 3 , C4 , C5 , C6  ,CZ,C8,C9,Clfl,Cll,C12 
DIMENSION  R(M1  ,M1)  ,S  (Ml,  Ml)  ,T  (Ml,  Ml)  ,G  (Ml,  NRHS) 

AL  = 1 . 

IF(LN.EQ.l) AL=0 . 

UO  1 11=1, K15 
G (It, NRHS) =0. 

00  1 Jl= 1 , K15 
R( II, J1)=0  . 

S(I1, J1)=0. 

T (I1,J1)=0. 

1  CONTINUE 
J1=0 

JO  2 Il=K8,K15 
T ( I 1 , Il)=l. 

Jl=Jl»l 

R(I1,J1)=-1. 

2  CONTINUE 

EQUILIJRIUM  EQUATION  FOR  1=0 
R(1  ,K8)  =D0*rUH-Qii**  2/011 
T (1,K6)=2.*911/  (RR*D11)  *XNXX 
T(1,1)=1./(RR**2*D11) 

G ( 1 ,NRHS)  = -XNXX*AWZPP(  JP  , 1 ) «■  AXPRE  S ( JP  , 1)/ 

E4  = -C4/2 « 

E5  = -C 5/2. 


E6=C2/2. 

00  3 J= 1 * KFQUR 
JS=J**2 


Jl  = J*1 

T ( 1 ,K3*J)  = T (1»K3*J)  +JS*(E4*(AL*AWM(JP,  Jl)  *AWZ( JP,  JIM 
lfrE6*Ai_*CFM(JP,J)  ) 

T ( 1 , J i ) = T ( 1 , J1 ) + JS  * (E4*( AL*AWMFP( JP,  Jl) f AWZPP(JP, Jl) ) + E5* (AL* 
1AWM(JP, Jl) *AWZ( JP, Jl)  > +E6*AL*CFMPP ( JP, J) ) 
S(1,J1)=S(1,J1)+2.*JS* (E4* (AL*AWMP(JP, Jl) * AWZP ( JP , J 1) ) «-E6*AL* 
1CFKPI JP, J) ) 

T (1,K9+J)=T (1,K9+J) + JSME4* (AL*BWM(JP, J ) + 8 WZ ( JP , J ) ) +E6 * AL* 


10FM ( JP, J) ) 

T ( 1 , K 1 + J)=T ( 1,  K 1+J) +JS*(E4*< AL*awMPPCJP,J) *BWZPP(JP,J)  ) *E5*(AL* 
13WM  (JP,  J)  4-OWZ(  JP,  J)  ) *E6*AL*0FMPP(JP,  J)  ) 

S(1,<1»J)  = S(1,K1*J)  + 2.*JS*  (E4*  ( AL*3WMP  (JP,  JH-BWZP  (JP,  J > )«■ 
1AL*£o*0FMP( JP, J ) ) 

T( 1,K11*J)=T(1, K11*J> +JS*E6* (AL*AWM( JP.Ji) +AWZC JP, Jl)) 

I ( 1,K13+J)=T (1 , K13  + J)=JS*E&MAL*3WM(  JP,J)+3WZ(  JP,J)  > 


T(l»K3*J)  = T(l»K3  + J)*JS*EbMAL4AWMPP(JP,JlJ *AWZPP ( JP , Jl )l 
T (1  , K 5 + J ) = I ( 1 , K 5 * J)  *JS*Eb*  (AL*UWMPF( JP, J)  fBWZPP  (JP,  J)> 
S(1,K3*J)=S(1,K3*J> ♦Z.*JS*£b» (AL* AWMP(JP, Jl) ♦ A W ZP ( JP, J 1 ) I 
S(l,K5*J)=S(l,K5*J)*2.*JS*E6MAL*dWMP(JP,J>*BWZP(JP,JI> 

G ( 1 ,NRHS>=&(  1,NRHS)  »AL»JS*  (E4*  l AWH  ( JP,  Jll  * AWI1PP  (JP,  Jl)  ♦ 
1AWMP(JP,J1)**2*BWM(JP,J) *BWMPP ( JP , J) fBWMP(JP.J) **2>+£5/2.* 

2(AWM( JP, Jl) **2«-DWM (JP, J) **2)  ♦ tb* ( ArfH ( JP , J 1) *CFMPP ( JP, J) ♦ 

3JWM (JP,  JDOFNPP  ( JP,  Jl  »2.  *AWMP(  JP,  Jl) *CFMP( JP,J) *2  # *BWMP ( JP , J I 
4*QF  MP  ( JP  , J)  *AWMPP(  JP,  Jl)  *CFM  (JP,  J)  +DWMPP(  JP,  Jl  *DFM(  JP,  Jl)  I 
3 CONTINUE 

EQUILIBRIUM  EQUATIONS  1=1,2,,,  KFOUR  SET  A 

DO  4 11  = 2,  K1 
1=11-1 
IS=I**2 
IS2=IS**2 

R ( 1 1 , I1K8-1)  = D0*H11 
R!I1,I1*K11-1)=-Q11 
T(Il,Il*K8-l)  = -Cb*IS«-XNXX 
T (11,  IDK11-1)  = C7*IS-1  ./RR 

T (II , 11) =C8*IS2-Al*IS* (C4*AWMPP  ( JP , 1 ) *C5* A WN ( JP , 1 ) I 
IF(LN.EQ.l)  GO  TO  51 
DO  5 J=l, KFOUR 
J1=J*1 

T ( II,  III  =T  (II,  III  i-AL*C9*IS*J**2M  (AMM(  JP,  Jl)  ♦2.*AWZ(JP  ,J1)  )* 
1AWM(JP,  Jl)  M8WM  ( JP,  J)  «-2.*6WZ  (JP,J)  )*OWM(JP,J)  I 
5 CONTINUE 
51  CONTINUE 

T ( 1 1 , I1«-K3-1)  = -C1C*IS2 

S(I1, I1*R1-1)=-2.*XNXY*C1*I 

T(  I1,K8)  =-C4*IS*  ( AL*AWM(JP  , ID  + AWZ1JP,  ID  ) 

Till,  1)=-C5*ISM  AL*AWM(JP,  II)  +AWZ  (JP.Il)  > 

G (II  , NRHS)  =-XNXX*AWZPP  ( JP,  ID  ♦ 2 . * XUX Y * Cl  * I *bKZP  ( JP  , I)  ♦ AXPRES  ( JP, 
111)  -C4*IS*AL*AWM  ( JP,  ID  *AWMPP  l JP,  1)  -C5*IS*AL*AWM(JP,  I1MAWM  (JP,1I 
E'5  = C9*IS*AWM(JP  , ID  *AL 
E4  = C9*ISMAL*AWM(JP,I1)*AWZ(JP,  III) 

DO  6 J=l, KFOUR 

J1=J*1 

JS=J**2 

IF (LN.EQ.l) 'GO  TO  7 

Gill, NRHS) =G (II ,NRHS) *E4*JS* l AW M ( JP , J 1) ** 2+&WM ( JP , J ) ** 2 ) +E5* JS* 
lt( AWM (JP  , Jl) *2. *AWZ( JP , J 1) ) *AWK (JP, Jl) +(BWM(JP, Jl *2.*3WZ(JP, J)l * 
23W  M ( J P , J ) ) 

7 T(I1,J1>  = T(I1,J1D2.*E4*JSMAL*AKM(JP,J1)  *AWZ(JP,Jl)l 
T(Il,KD-J>=T(Il,Kl+J)»2.#E4*JS*(AL’BHM(JP,J)«-BWZ(JP,j)) 

6 CONTINUE 
E5=C2/2. 

00  3 J= 1 , KK2 

Till, Kll »J) = T(I1, K1DJ)  + E6M  ALB  (I  ,J,1, AW Z , JP , N J , N AW , 1 , 1 ) 
1*’ALB(I,J,4,AWZ,JP,NJ,NAW,1,1)I 
T(Il,K13*J)=T(Il,K13*J)+E6MALa(I,J,l,BWZ,  JP,  N J , NBW , 1 , 2) 
l»ALd(I, J.7,3WZ, JP,NJ,NBW,1,2)I 
S(I1*K3  + J)=S(I1»K3  + J)+£6MAL«3(I»J»3,AWZP,JP»NJ»NAW,1,1) 

1» ALB( I , J , 6 , AWZP , JP, N J , NAW , 1 , 1) ) 


S(Il»Kb*J) = S (II ,KS*J)  + E6* t ALUti , J , 3 , dW 2P, JP, NJ , NOW , it  2 1 
1*AL3< I ,J,9, 3WZP.JP. NJ.NBW, 1,2) ) 

I <Il,K5*J)=TlIl,<i*J)*EbMAL8tI  , J ,2 , AW  29? , JP, N J , NAW , 1 , 1 1 
1*ALG< I, J,5,AWZPP, JP.NJ.NAW, 1,1)  I 
MIl,K5*J>=MIi,<&*JI*E6MALBtI  , J,2,  BWZPP,  JP.NJ  ,N6W,  1,  2) 

1 ♦ A L 3 f I,  J, 3, BWZPP,  JP,NJ,NBK,1  ,2)  ) 

IFtLN.EQ.l)  GO  TO  8 

T ( Il,Kil*J)=T  ( II  , K1 1 * J ) *E6*  (ALB  ( I » J , 1,AWM,JP,NJ»NAW,1»  1 ) * 
lALdtl ,J, 4, AWM, JP.NJ.NAW, 1, 1) > 

TlU,K13*J>  = TlIl,K13*J)*CGMALBt  1 , J . 1 , BWM  , JP , N J , N8W  • 1,  2)* 
lALJtl  ,J,  7.BWM,  JP,NJ,  NBW,  1, 2)  ) 

St  II,  K3*  J)=StIl  ,K3*J)*E6*tALBlI,J, 3,  AW  MR, JP.NJ.NAW, 1,1)* 
1ALB(I,J,6. AWMP, JP.NJ.NAW, 1,1) ) 

S t II, K5*J) =StIl ,K5*J) *E6* ( ALBt I , J,3,dWMP, JP,NJ,NBW, 1,2) ♦ 
lALb(I,J,9.BWrtP, JP,NJ,NBW,1,2) » 

T tIl,K5*J)=T  til  » K3  * J ) *E6MALBtI  ,J»2,AWMPP»JP,NJ,NAW»1,  1)  * 

1AL3II , J,  5,  AWMPP,  JP.NJ.NAW,  1,1)  ) 

TtIl,K5*J)=TtIl,K5*J)*E6MALB(I  , J , 2 , UW  MPP , JP,  N J ,NBW  , 1,  2)* 
lALBtl, J, 8, BWMPP.JP, NJ.N3W, 1,2) ) 

IJ1=I*J 

IF1IJ1.GT.KFOUP)  GO  1 0 8 0 

T(I1,IJ1*1)=T1I1,IJ1  + 1)*E6*'(ALB(I,J,1,CFMPP,JP,NJ,NF,2,2)I 
T<  II,  IJ1*K1)=T  t II,  IJ1*K1)  *E6*AL3t  I,  J,  1,  OF  MPP,  JP.NJ.NF,  2,2) 

Stll,  Ulfl)=StIl,IJl  + l)«-E6*ALe<  I,  J,3,CFMP,  JP,NJ,NF,2,2I 
St  II, IJ1+K1)=S  t II, IJ1+K1) +£6*AL3t I,J , 3,0FMP, JP, NJ,NF,2 ,2) 

Tt  I1,IJ1*K8I  =T  t Il.XJl+Kft)  *Eb*AL3t  I,  J , 2 , CFM , JP,  N J,  NF, 2,  2) 

Mil,  IJ1*K9)=T ( II, IJl*K9)+E6*AL3tI,J,2,DFM,jP,NJ,NF,2,  2) 

80  IJ2-IAUS (I-J) 

IFt  IJ2.GT. KFOUR)  GO  TO  90 

Till,  IJ2H)=TtIl,lJ2*l)*E6*ALBt  I , J , 4 , C FMPP , JP , N J,  NF, 2,  21 
SlIl,IJ2*l)=SlIl,IJ2*l>*Et>*ALBt  I , J , & , C FMP , JP , N J , NF  , 2 , 2 ) 
ft  II,  IJ2*K3)=T  1 11,1  J2  + K3  )«■£&*  A LJ<  I,  J,  5,  CFM,  JP.NJ.NF,  2,2) 

IF  ( U2.EQ.  0)  GO  TO  90  , 

MI1,IJ2*K1)=T(  11,1  J2*Kl)*E6*ALOt  I,  J,  7,  OFMPP.JP,  NJ.NF,  2,2) 

St  II, IJ2*K1>  =S t II, IJ2+K1) +£6*AL3t I, J ,9, OF, IP, JP,  NJ.NF, 2 ,2) 

Mil,  IJ2+K9I  =T  (11,1  J2+K9) *E6*AuBt  I,  J.tt.DFM,  JP,  NJ,NF,2,2I 
9 0 Glll,NRHS)=G(Il,NRHS>*E6M  tALBt  I,  J,l,  AWM,  JP*NJ,NAW,  1, 1)* 

1AL it  I ,J, 4, AWM, J?,NJ,NAW,i, 1) ) *C FMPP t J P , J) * t ALB  1 1 , J , 1 , B WM,  JP , N J, N3W 
2,1,2)  *ALBII» J,  7, BWM.JP, NJ.NBW, 1,2))*  OF  MPP  t JP , J ) *(ALBtI,J,3,AWMP 
3, JP.NJ.NAW, 1,1)  *AL8tI»J»6»AWMP, JP.NJ, NAW, l,l))*CFMPtJ?»J)* 

4 1 AL  J t I,  J,3,BWHP,  JP,  NJ.NBW,  1,2)  + AL6(I  , J , 9,  BWMP  , JP , NJ  , No  W , 1 , 2 ) > ♦ 
BOFMPt JP, J) *t AL 311, J, 2, AW MPP, JP.NJ.NAW, 1,1) + ALB l I , J , 5 , A WMPP , J P, N J 
6,  NAW,  1 , 1))*CFH  ( JP,  J)  * t ALBt  I,  J,2  ,3  WMPP,  JP,  N J,  N3W  , 1 , 2) 

7tAt_3(i,  J,  8,  BWMPP.JP,  NJ.NBW,  1,2)  )*OFM(JP,J)  ) 

& CONTINUE 
4 CONTINUE 

C EQUILIBRIUM  EQUATIONS  1=1,2,,,,  KFOUk  SET-  d 
1=0 

OO  14  I1=KZ,K3 
1=1*1 
111=1*1 
IS=  I*  *2 
I32=IS**2 

I 


1 


R{  1 1 » K9 ♦ 1 ) =l)J*  H 1 1 
RII1,K13*II=-Q11 
rill,  K9*I)=-Cfc*  IS*XNXX 
T l 1 1 »K13*I )=C7* IS-l./RR 
I ( 1 1 * <5 ♦ I ) --C1 0 *IS 2 
Sill, I11)=2.*XNXY*C1*I 

T ( II ,Kb) =-C4*lS* IAL*BWM( JP,I) + 8WZ (JP,I )) 

Till, 1)  = -C5*IS* ( AL'dWM  CJP, I) + BWZ (JP, I) ) 

Till, 1 ♦ K 1 J =Cd*I S2-C4*IS*AL*AWt',PP( JP, 1) -C5*  IS*AL*AWM (JP , 1) 

G ( Il.NRHS)  =-XNXX*BWZPP  (JP,  I)  - 2.  *XNXY*C1*I*AWZP  (JP,  Iill  * 

1 3XPR- S ( JP, I) -C4*IS*AL*BWM(JP.I> *AWMPP( JP* 1 I -C5* IS* AL*B WH l JP, I) 
2AWM  (JP,  1) 

E4=C9*IS*( AL*BWM( JP, I) *BWZlJP,III 
£5  = C9*IS*AL*t?WM ( JP, I) 

00  15  J=1 , KF OUR 
JS=J**2 

ji=j-h 

IFILN.EQ.ll  GO  TO  17 

T (I1,K1*I)-T(I1»K1*I) *C9*IS*JS* l (AWMIJP,  J1)*2.*AWZ(  JP,  Jl)l  * 
1ANM  ( JP,  J 1)  KBWjI  ( JP,  J)  +2.  *BWZ  (JP  , J)  )*BWM(JP,J)  ) 

Gill,  i'lkHS)-G(Il  , NRHS)  *E4*JS*  ( A W M ( JP , J1 ) **  2 + BWM  ( JP,  J)**  2)  + 
lE5*JSM(AWM(jP,Ji)*2.*AWZ(JP,Jl))*AWMlJP,Jll*(3WMlJP»JI*2.* 

2GK  L ( J P , J ) ) * J AM  ( J P , J ) l 

17  Till, J1)=T(I1, Jil*2. *E4*JS*lAL*AWM(JP, Jll  *AWZ(JP, Jill 
T ( 11, K1 + J) =T (II , <1  ♦ J) +2.*E4*JS* l A L* BUM (JP,J)*3WZ(JP»J) I 

15  CONTINUE 

E6-=C2/2  • • 

00  18  J= 1 , KK2 

THl,Kii*J)  = T(I l,Kll*J)*E6*(ALB(I,J,l,eWZ,JP,NJ,N3W,l,2>- 
1AL31I ,J, 7.BWZ, JP,NJ,N3W, 1, 2) I 
Till, K1 3 +J)=T( I 1,K13*JI + E6M-AL5II, J, i,AWZ ,JP,  NJ,NAW,1 , 1) ♦ 
1ALBII , J.4.AWZ,  JP.NJ.NAW.l, 1)1 
S(  I1,K3*J)=SII1  ,K3*J>  + E6*  (ALB (I  , J, 3,  BWZP, JP.NJ ,NBW, 1,21- 
1ALBII ,J, 9, BWZP.JP, NJ,MBW, 1,2)) 

S(I1,K5+J1=SII1,<5*J) +E6*(-ALB(I, J.3.AWZP, JP,NJ,NAW,1, 11* 
1AL01I ,J,6, AWZP,  JP.NJ, NAW, 1,1) ) 
T(I1,K2*J)=T(I1,K3*J)*E6MALBII,J,2, BW  ZPP  , JP,  N J , N 3W  , 1 , 2 1 - 
lALOd  , J,  8, 9WZPP  , JP,NJ,N3W,1,2>) 

T (II  ,K5*J)  =T  (II  , K5+J)  + E6M-AL61I,  J,2,  A WZPP , JP , N J , NA W , 1 , 1)  * 

1 AL8(I,J,5,AWZPP,JP,NJ,NAW,1,1I) 

IF(LN.EQ.l)  GO  TO  18 

T ( 1 1, <1 1 +J ) = T ( I l,Kll*JI*E6*l ALB (I,J, 1, BWM , JP, NJ , N9W , 1, 2J- 
1ALJII , J,7,3WM, JP,NJ,NBW, 1,2)  1 
T III  ,M3*J)  = TII1,K13*J)+E6*  I-ALBII,  J,1,AWM,JP,  NJ.  NAW.l  , II* 
1ALQII,J,4,AWM,  JP.NJ,  NAW,  1,11  1 
S(I1,K3*J)=S(I1,K3*J) *E6* ( ALBl 1 , J , 3, BWMP , JP.NJ , NBW , 1, 21 - 
1AL31I  , J,9,BH'1P,  JP,NJ  , N BW  , 1 , 2 ) ) 

S(I1,K5*J)=S(I1  ,.<5*J)  *E6* l-ALB (I , J,3, AWMP, JP.NJ  ,NAW*1, 11* 

1AL91 I ,J,6,AWMP,  JP,NJ, NAM, 1,1)1 
TII1,K3*J>=T(I1,K3*J)*E6*IALU(I,J,2,8WMPP,JP,NJ,N3W,1,21- 
1ALSII , J, 3, BWKPP  , JP.NJ, NBW, 1,2)) 

T( II ,K5*J)-T III ,K5  + J) *E6* l-ALB  1 1, J, 2, AWMPP , JP, N J, NAW, 1 , 11 ♦ 

1 ALb (I, J,5» AWMPP, JP.NJ, NAW, 1,1)  I 


r ^ 


IJ1=I*J 

. IF  (IJ1.GT.K.F0UR)  GO  TO  81 

T Ml.  IJlfrKl)=T(  II,  IJ1+K1)  «-E6*ALlJ(  I,  J,1,CFHPP,  JP,NJ,NF,2,2) 
H11,IJ1*1)*I(  1 1.  IJl  + l)-C6*ALO(  I,  J,l,OFMPP,JP,NJ,NF,2,  21 
. S(Ii,  IJ1«-K1)=S(  I1,IJ1+K1)*E6*ALB(I,J,  3,CFHP,  JP,  NJ,NF,2,2) 

S (II,  IJ1+1)=S(I  1 ,IJlH)-Eb*ALB(  I,  J,3,DFMP,  JP,  N J , NF  , 2 » 2 ) 

Till, IJ1+K9)=T( I1,IJ1+K9)+E6*ALB(1,J,2.CFM,JP,NJ,NF,2,2) 

T <Il*IJl+K8)=T ( Ilf  I J1+  K8) -E6*AL8( I , J , 2 , OF M , JP, N J , NF , 2 , 2 ) 

81  IJ2=IA0S  (I-Ji 

IF( IJ2.GT. KFOUK)  GO  TO  91 

Mil,  IJ2  + i)  = nil,IJ2  + l)  + E6*AL6(I,J,4,0FMPP,JP,NJ,NF,2,2> 

SCI1,  IJ2+1)=S(  1 1,IJ2*1)+E6*ALB(  I,  J,6  , OFMP,  JP,NJ  ,NF  ,2,2) 
T(Il,IJ2+K8)=T  (U,IJ2  + K8)+E6*ALB(  I,  J,5,DFM,JP,  NJ,NF,2,2) 

IF ( I J 2 • EQ • 0 ) GO  TO  91  ' 

T ( II, IJ2+K1)=T ( I1,IJ2+K1)-E6*ALB( I, J , 7 , CFM PP, JP , N J , NF, 2,21 
S(I1,  IJ2*K1)=S(  II, IJ2+K1)-E6*ALB( I, J.9.CFMP, JP, NJ,NF,2 ,2)  ' 

T(I1,IJ2«-K9)=T<I1,IJ2+K9)-E6*Al8(I,J,8,CFM,JP,NJ,NF,2,2* 

91  G ( II,  N?.HS)=G  (II  , NRHS)  + E6*  ( (AL9  ( I , J , 1 , BWM , JP, N J , NOW , 1 , 2 ) - 
1AL  B 1 1 ,J, 7, BWM, JP,NJ,NBW, 1, 2) ) * C FMPP ( JP , J) 

2M- AL  3(1, J,1,AWM,JP,NJ,NAW, 1,1)  +ALB < I , J , 4 , AWM , JP , N J , NA W, 1 , 1 ) )*DFM 
4PP ( JP , J ) *■  ( ALb(I  ,J,3,BWMP,JP,NJ,  N3W,1,2)-AL3(I,J  ,9»OWMP  ,JP,NJ,NUH, 
51,2) ) *CFKP  ( JP,  J)  +1-AL8  (I  , J , 3 , A W MP  , JP  , N J,  N A W , 1 , 1 ) + ALB(  I , J , 6 , A WMP ,, 
6JP ,NJ,NAK,1,1) ) *DFNP( JP, J) MALB (I  , J , 2 , BWM  PP,  J?  , NJ  , NOW  , 1,2)- 
7AL3(I  ,J,d,3WilPP,  JP,NJ,  NBM,  1, 2) ) *CFM ( JP  , J)  + (- ALB  (I,  J,2,  AWHPP, 
8JP.NJ,  NAW.l,  1)  <-ALQ(I,J,5,A  KMPP  , JP»NJ,NAW»  1 , 1 ) ) * OF  M ( JP,  J ) ) 

18  CONTINUE 
14  CONTINUE 

C COMPATIBILITY  EQUATIONS  1= 1 , 2 , , , 2 *KFOUR 
1-0 

00  19  I 1=K4 ,K5 
1*1*1 
IS=I**2 

R(I1»K11+I)=011 
T(I1,K11+I)=-C11*IS 
T(il,K3+I)=C12*IS**2 
IF (I. GT . KFOUR)  GO  TO  35 
R(I1,  K8H)=Q11 
T ( 1 1 » K8  + I)  = -C7*  IS+l./RR 

T(Il,I+l)=Clu*IS*+2-C2/2.*IS*(AL*ANMPP(JP,l)) 

T II1,K8)=-C2/2.*IS» (AL*AWM(JP,I+1)+2.*AWZ(JP, 1+1)1 
G ( I 1, NRHS) =-C2/2 • * I S* ( AL* A WM (JP , I + 1) * A WMPP ( JP ,1)1 
35  E6--C2/2. 

E5=E6/2. 

DO  2D  J = 1 » KFOUR 
Jl=J+i 

T (II , K8+J)  -T  (II  , K8+J)  + E6MALBI  I , J,l,  AWZ,  JP.NJ,  NAN  , 1,1)  ♦ 

1AL  J ( I »J»4,AWZ,JP,NJ»NAW,1,1)  ) 

T ( II, K9+ J ) - T (II , K9  + J ) + E6* ( AL d ( I , J ,1,  BWZ,JP,NJ,  N3W,  1,2)  ♦ 
lALBd  , J,  7,  BWZ,  JP,NJ,N8W,  1,2)  > 

S(I1,  J+1)=S(I1,  J + l) +E6M  ALBli,  J , 3 , AW  ZP  , JP  , NJ,N  AW,  1 , 1) 

1 + A L b ( I, J.6.AWZP, JP,NJ,NAW,1,1) ) 

S ( I1.K1+ J)*S (11 ,K1+J) +E6* (ALB(I, J,3,  3H ZP, JP, NJ , NSW , 1 , 2 > ♦ 
lALod ,J,9,BWZP, JP,NJ,NBW,1,2)) 


r -1 


I (1 1,  J*1)  = T ( II,  J*l)  + E6M  ALb(  I,  J ,2  .AWZPP,  JP.NJ,  NAW  ,i  ,1)  «■ 

1ALB  (1 , J, 5. AWZPP, JP, NJ, NA W,i, 1) ) 

T(il,Kl*J)=MIl  ,K1+J)  »E6*  (ALB  (I  , J ,2,3WZPP,  JP,NJ,N3W,1,  2)* 

1ALBU  ,J.  8,  BWZPP  , JP.NJ,  NBW,  1,2)  ) 

if  tLN.ta.n  go  to  20 

T(Il,K8*J)  = T(Ii,K8*JI*E5MAL8(I,J,l,AWM,dP,NJ,NAK,l,l)«- 
IAlD  ( I .Jf  4,  AWN,  JP,NJ,NAW,  1, 1)  ) 

I ( il,K9*J)  =T  (II  ,K9*JI  *£5MALB<  I,J,1,BWM,  jf»  ,NJt  NBW  » It  2)  ♦ 
1AIB<I,J,7,BWM,JP,NJ,NBW,1,2‘>) 

S(Il,  J*1)=S(I1,  J + 1)*E5MALB(I,  J,3,AWNP,JP,NJ,NAW,1,1>* 
lALbd  tJt  6,  AWNP,  JP.NJtNAWtltl)  ) 

SI  I1,K1*J)=S(I1  ,<1*J)  *E5MAL3(  I , J,3,6WMP,  JP,NJ,NBW,1,2»  ♦ 

1ALBII ,J,9, EWMP,  JP,NJ,NBW,1,2>) 
rill, J+1)=T(I1, Jf 1) +E5*( ALBI I, J,2,AWMPP,JP,NJ,NAW,1,1)  ♦ 

1 At  3 ( I ,J, 5, AWMPP, JP,NJ,NAW,1, 1)  ) 

T (I1,K1*J)=T  (I1,K.1  + J)*E5*(ALB(I  , J , 2,  BW MPP , JP , N J , NBW , 1,  2>fr 
1AL3II,  Jt  3,  BWKPP.JP,  NJ.NBW,  1,2)) 

IJ1=I»J 

IFIIJI.GT.KFOUR)  GO  TO  36 

T ( lit  IJ1H)  = TI  1 1 ,IJ1  + 1)*E5*ALB(  I , J , 1 , AWMPP , JP , N J,  NAW,2  , 1» 

T (11,  IJ1*K1)=T  < Il,IJl  + Kl)+E5*Ai.3d,  J,1 ,0WMPP,  JP.NJ.N3W  ,2,2) 

Sill,  IJ1+1)=S(I 1,1  JH-1)*E5*AL6  ( I, J,3, AWHP, JP.NJ.NAW, 2,  1) 

SI  II, IJ1*K1)=S( II,  IJ1  + K1) *E5*AL3I I, J , 3 , BWMP, JP . NJ , NBW, 2,2) 
T(Il,IJl*Ktt)=Tt Il,IJl*K8)*L5*ALBlI,J,2,AWM,JP,NJ,NAW,2,i> 

I ( II , IJ1»-K9>=T ( II, IJ1*K9) *E5*AL3( I, J. 2,3WM,JP,  NJ,NBW,2,2) 

36  IJ2=I A8S  C I — J ) 

IF  ( IJ2.GT.KF0UR)  GO  TO  37 

Till, IJ2*l)  = nil,IJ2*l)*E&*ALBII,J, 4, AWMPP, JP.NJ.NAW, 2,1) 
SII1»IJ2*1)=SII1»IJ2*1)*-E5*AL&II,J,6,AWMP,JP»NJ,NAW»2»  1) 
Tdl»IJ2*-K6)=T(Il»IJ2*i<8)*E5*ALB(I,J»5»AWM»  JP.NJ.NAW, 2,1) 

IF ( I J2 .EG. u ) GO  TO  37 

Ttll,  IJ2*K1)=T  ( IlfIJ2  + Kl)*E5*ALBlI,J  ,7,  3WMPP  » JP , NJ  » NBW  ,2,2) 

S(  II,  I J2  *K1 ) -S  ( II,  IJ2  + K1)  *£5*AL3(I,J  »9,-3WMP,JP»  NJ  » NBW  , 2,2) 

Till, U2  *K9)  =T  III, IJ2*  K9  )*E5  * ALB(I,J, 8, 8WM,JP,NJ, NBW, 2, 2) 

3 7 Gl II, NRHS)  =G  (1 1 , NRHS)  *E5* ( (ALBl I,J,1,AWM,JP,NJ,NAW,1,1)+ 

1ALBII ,J, 4, AWM, JP.NJ.NAW, 1,1) )*AWMPP(JP,J1>  M ALB (I ,J,i,  BWM,  JP , N J , 
2NSW, 1 ,2) *AL3(I,  J ,7,OWM, JP,NJ ,NBW, 1,2) ) *BWMPP(JP, J) ♦ ( ALB(I,  J,  3, 
3AW  3P,  JP.NJ.NAW,  1, 1M-AL3I  I,  J.6,  AWMP,  JP.NJ.NAW.l,  II  I *AWMP(JP,  Jl)  ♦ 

3 ( Al J ( I , J , 

43,6WMP,JP,NJ,NBW,1,2)  f AL B ( I , J , 9 , 3WMP , JP , N J , NBW , 1 , 2) ) *3 WMP ( JP  , J ) 
5(ALB(I,J, 2, AWMPP, JP.NJ.NAW, 1,1)  *AL3(I, J,5 , AWMPP  ,JP,  NJ  , N AW  , 1 , 1 ) ) 

- 6AWM  (JP,  Jl)  MALB  (I,  J,  2,  BWMPP.JP,  NJ.N3W,  1,2)  *AL9  ( I,  J . 8 , 8 WMPP  , JP  , 

7N J, Nd W,1 ,2) ) *3KM  (JP, J)  ) 

20  CONTINUE 
. 19  CONTINUE 

C COMPATIBILITY  EQUATIONS  1=1,2,,,,  ,,2*KFOUR  SET  B 
1=0 

00  21  Il=K6,K7  ' 

1=1*1 

IS=I**2 

R(I1,K13*I)=D11 
T (I1,K13*I) =-Cl 1*IS 
T(I1,K5*I)=C12MS**2 


“Tl 


im 


IK  ( I.  GT  . KFOUtfl  GO  FO  45 

«mi,<9«-i)=oii 

T-(  II  , K9*I)-~C7*  IS*-1./RR 

T(U,Kl*-I)-Glu*IS**2-C2/2.*IS*(AL*AHMPP(JP,l)l 
T(  II,  n(J»=-G2/2.  ♦ISMAL*3WM(JP,I)«-2.*3WZ(JP,I)» 

’ G ( Il,Ni\HS)  =-C2/2.*lS*  ( AL*bWU  (oP  , I ) *A  WMPP  ( JP,  1)  1 

45  E6=-C2/2. 

E5=E6/2. 

00  22  J= 1 » KF OUR . 

J1  = J*  1 

T(  I1.K3  + J)  =T  (II  ,K8*J)  *E6MALB(  I ,J,1,  BWZ,JP,NJ,  N3W,1,2)  - 
1AL  «(I,J, 7.BWZ, JP.NJ, NOW, 1,2) ) 

T ( 11, K9*J)=T (II , K9*J)  *E6M-ALG(  I , J,1  , AWZ,  JP.NJ  ,NAW,l,t  )♦ 
1ALJCI,J,4, AWZ.JP, NJ.NAW, 1,1) ) 

S(I1,J1)=S<I1,J1)  *£6* ( ALB( I, J, 3 , 8WZP , JP, NJ , N8W , 1, 2 ) - 
1AL'3(I  ,J,9,  BWZP,  JP,NJ,NBW,1,2)  ) 

S(I1,K1*J)=S(I1 » K1 ♦ J ) + E6* (-ALB( I, J, J, AWZP, JP,NJ,NAW,i, 1)* 
1ALJII ,J,6,AWZP,  JP,NJ,NAW,1,1)) 

T ( II,  J1)=T(I1,  J1)+E6MAL3(1,J,2,8WZPP,  JP,NJ,N3W,1,2)- 
1AL3  (I , J,  S.tWZPP  , JP,  NJ,  NOW,  1, 2)  ) 

T ( II ,K1*J) =T (II ,K1*J) *E6» (-ALB( I, J,2,AWZPP,JP,NJ,NAW,1 ,1)+ 
lALud  ,J, 5,  AMZPP.JP,  NJ.NAW,  1,1)  ) 

IF(LN.EQ.l)  GO  TO  22 

T<I1,K8*J)=T (Ii,<8*J)  «-E5MALB(I,J,l,8WM,  JP.NJ,  NOW,  1,2)- 
1AL»3(I  ,J,  7.8WM,  JP.NJ.NOW,  1,2)  ) 

T(I1,K9«-J)-T(I1,K9  + J)«-E5*(-ALB(I,J,1,AWM,JP,NJ,NAW,1,1>* 
1AL31I, J, 4, AWM.JP, NJ.NAW, 1,1)) 

S(I1,J1)=S(I1,  J1)*E5*  (ALL>(I,J»3,dWMP»JP,NJ»NBW,l,2)  - 
lALOd  , J,  9,bWMP,  JP,NJ,NeW,l,2)  ) 

S ( 1 1 , K 1 ♦ J)  =S  ( 1 1 ,K1  + J)  *E5*  (-ALB  (I,  J,3,  AWMP,  JP,NJ,NAW,1,  1)  *• 
lALBd  , J,6,AWM°,  JP,  NJ.NAW,  1,1)) 

T ( 1 1 , J1 1 = T (Il»Jl)*E5* (ALB(I,J,2,BWMPP,JP,NJ,N3W,1,2)- 
lALd(I,J,8,bWMPP,  Jt*,NJ,N3W,l,2)  ) 

T (I1»K1+J)-T  (Il.Kl+J)  ♦ E5M-ALB(  I , J , 2 , A WMPP  , JP,  N J , N A W , 1 , 1)  ♦ 
lALa  ( I , J,5,AWrlPP  , JP,  NJ.NAW,  1,  1)  ) 

IJ1=I*J 

IF (IjI.GT.KFOUR)  GO  TO  46 

T(I1,IJ1*K1>=T(  II, I J1  + K1)*E5*ALB( I, J ,1,AWMPP, JP , NJ.NAW ,2, 1) 
T(I1*IJ1*1)=T(11,  IJ1+1)-E5*ALB(  I,  J , 1 , B WMPP  , JP  , N J , N BW  , 2 , 2) 
S(I1»IJ1*K1)=S  (II, IJ1*-K1)+E5*AL9(I»J»3»AWMP,JP, NJ.NAW, 2,1) 
3(11,  IJ1*1)=S(  11,1  JH-D-Eb* ALB  ( I, J,J,BWMP, JP , N J , N 8W  , 2 , 2 ) 
T(I1,IJ1*K9)=T  ( I1.IJ1+K9)  *E5*AL3(  I,J,2,  AWN  , JP  , N J,  NAW  ,2 , 1) 
T(I1,IJ1*K3)=T( II,  IJ1+K8) -E5*ALB( I, J,2,BWM,JP, NJ,N8W,2,2) 

46  IJ2=1 AdS (I-J) 

IF(  IJ2.GT.KF0UR)  GO  TO  47 

r(Il,IJ2*l)=T(Il,IJ2*l)*E5*AL8 ( I , J, 4 , B WMPP , JP  , N J , NOW , 2 , 2) 
S(I1,IJ2H)=S(1 1,IJ2*1)+E5*ALB(I,J»6»BWMP,JP,NJ  ,N0W,2,  2) 
T(Il,IJ2»-K8)=T (11,1 J2+K8)*E5*ALB(1, J , 5 , BW M , JF , N J , N BW , 2 ,2) 
IF(IJ2.£Q.0)  GO  TO-  47 

Till,  IJ2*K1)=T  ( 1 1,IJ2+K1)-£5*ALB(  I , J , 7 , AW  MPP  ,JP,  NJ.NAW  ,2, 1) 
SI  II , IJ2*K1) =S ( II, IJ2+K1)-£5*ALB(I,J , 9 , AW  HP, JP , NJ , NAW , 2,1) 

T (II, IJ2  + K9) = T ( II, IJ2+K9)-E5*AL3( I, J,8, AWM, JP,  NJ.NAW, 2 ,1) 

47  G ( II,  NkHS)=G  (II,  NrfHS)  + E5M  ( AL8  ( I , J , 1 , 3WM,  JP,  N J , NBW  , 1 , 2 ) - 


72 


1 AL  t)  ( I , J,  7,uWM,  JP,NJ,NBW,  1,2)  )*  ANMPP(  JP«  Jll  ♦<-ALB(I»  J»1  . AWM.  JPt 
2NJ.NAW,  1,1)*  AL3  (I,  J,«o  AWM, JP,N J ,NAW, 1,  1)  ) * OWMPP  < JP,  J>  «■  C AL  B ( I , J, 
33*DWitP,JP»NJ,N3W«l,2)  - ALb<  I,J»'J,t3HMP,jP,NJ,NBW,l,2)  )*AWMP(JP,J1)*- 
4(-AL  > 1 1 , J,  3,  AWNP,  JF.fcJ  ,NAW,1 ,1M-ALB(I,  J,b.  AWMP  , JP,  N J , N A W , 1,1))* 

5 JWrtP  l JP,  J)  MALI)  Cl,  J,  2,  BWMPP.JP,  NJ.NJW,  1,2)  -AL3  t I , J ♦ 8 , B WMPP  , JP  , 

6NJ  ,NJW,  1 ,2)  ) *AWM  UP,  Jl)  M -ALB  11  , J,  2,  A WMPP  , JP , N J * N A W , 1 , 1 ) «■ 

7 ALt-f  I , J, 5,  AWHPP, JP,  NJ,  NAW,  1, 1)  ) *0W,1(  JP,J)  ) 

22  CONTINUE 
21  CONTINUE 
RETURN 
ENO 


SUBROUTINE  POTSN (POT.STRYU.STRAU, SI  RYG ,ST R AG , X NXX , XNXY I 
C STRYU  UNIT  END  SHORTENING  U FOR  Y = Q. 

CSTRAU  'AVERAGE  UNIT  ENJ  SHORTENING  U 
C STRYG  UNIT  ;NU  SHORTENING  GAMA  FOR  Y = C. 

CSTRAG  A Vi_  RAGE  END  SHORTENING  GAMA 

COriMCN/FOURl/KFOUR.Ri*  K2.K3,  K4 , K5 ,K6, K7.K8 ,K9 
COMMON/ FOUR 2/K10»KlltK12»K13»K14»K15  , K.K2 
COMM'JN/GEOM/RR,  DO.Hll,  H12,H22,U11,Q12,Q22,  Oil,  U12, 022 
C0,MMuN/FACT3R/C  1,C2,C3,C4,C5  ,Go,C7,C3,C9,ClO,Cll,C12 
COMMON/FACT2/0L 1 , D12  , DL J » DL 4 . 0 A 1 , D A2 , 0 A3 , 0 A4 , 0 B 2, DB3,DB4,XNI,tXXP 
CCMMUN/CINTG/NEOPOT.MI (500) 

COMMON/ F10FR/DELT A, AL1,GA1,AL2 , 3T2.GA2 
CONMON/FACT 3/XL , XH 

CCMMON/PRESl/AWM (100 ,3) , AWMP (100, 3) , AWMPP ( 100. 3) 

GO MM ON/PRES 2/ DWM (100*21 • RUMP ( 1 0 0 , 2) , BW MPP (100,2) 

COMM JN/PRES3/CF  M (1GO»4),CFMP(10G,4),CFMPP(100»4) 

COMMON/PRESA/OFM (100 , A) , DFMP (100, A) .DFMPP (100, A) 

GOMMON/PRES5/AW  Z (100,3) , AW ZP ( 1 0 0 , 3) , AW ZPP ( 10 0 , 3 ) 

COITMGN/PRESS/BWZ  (10  0,2),  3W ZP, < 1 G 0 , 2),  9WZPP(10B,2) 

_ COMMON/ X XL O A 0/ A XPRE S(100,3)  ,OXPRES(1GO,2) 

POT  = C . 

S T R Y U = 0 • 

ST  RAU=  0 « 

STRYG-0. 

STRAG=0. 

OU  10  11  = 1 « NfcQPOT 
E7  = 1 . 

IF  < II.  EQ.l.  OR.  II.  Eli.  NEOPOT)  £7  = 0.5 
El=-Ull*AWNPP(Il,l)-AWM(Il,l)/RR 
£2  = 0. 

00  11  J= 1 , KFOUR 
JS  = J*  *2 
J1=JH 

E2=E2+JS*  ( (AWM t II, Jl) +2.*AWZ(I1 , Jl) ) ♦AWHtll, Jll  +(8WM(I1,J)* 
12.*eWZ(Il,J) )*BWM(I1, J)) 

11  CONTINUE 

E 1=£ 1 +C2/A • *E2 

PE  1 = 1 ./Dli*tl**2*0Ll*AWMPMIl,l)**2-2.*XNXX*DA2/Dll*El-2.*AWM(Ilt 
11)*AXPkES(I1,1) 

E2=-OA2/Dll*£l-DA3*AWMPP(Il,i> 

PSYU=E2 

PSAU=E2*Q.5*AWMP(I1,1)*(AWMP(I1  , 1 ) +2 . * AWZP  (1 1 , 1 > 1 
“£3=0. 

DO  13  J1=1,K1 
OO  18  J 2 = 1 , <1 

l1=E1»AWMP(I1, Jl)*(AWMP(Ii,J2l+2.*AWZP(Ii, J2)) 

IF  ( Jl'.EQ.Kl)  GO  TO  19 
E2=E2*Jl*Cl*3HM(Il,Ji)*AWZP(Ii»J2) 

19  IF ( J2.EQ.K1I  GO  TO  18 

E3=E3  ♦ AWMP  ( 1 1,  J 1)  *J2*C1*  (bWM  (II , J2)  «■  BW  Z (II  ,J2)  ) 

18  CONTINUE 


74 


PSYU=PbY U+E 1/2. 

PSYG  = E2«-£3 
El  = 0 . 

£2=9. 

E3=  0 • 

£4=0. 

£5  = 0. 

£6=0. 

£3  = 0. 

£9=0  • 

£10  = 0 . 

30  2G  J = 1 , KF  OUR 
J1=J*1 
JS=J+* 2 
>JS2=JS4*2 

E1=£1*JS2*  (AWM(  Il»Jl)**2*bWH(Il  , J ) **  2) 

E2=£2*AWMPP(I1, Jl) **2*BWMPF(I1, J) **2 

£3=E3 »JS*  ( A Wild  1 ,J1)  *AWHPR(I1,  Jl)  4-BWM  ( II , J ) *BWM PP { 1 1 , J > > 

E4=E4  *JS*  (AHnPT  II,  Jl)  **2+CWMP(Il,  J)  **Z) 

£5=£5+AXPRES (11 i Jl)*AhH(IliJl) +6XPRES ( I1,J)*BWM ( 1 1 , J) 
£6=E6«-JS*AWM(I1,  Jl) 

E3=£3+AWMPPCI1, Jl) 

E9=£9«-AWMP(I1,  J1IMAWMP(U,J1)  + 2.*AWZP(I1,  J1))+BWHP(I1  ,J)M 
13W;1P(I1,  J)  *2.*bWZPtIl,  J)  ) 

E10=E10+JMAWMP  (II,  Jl)  *(3WM(I1,  J)  *BWZ(I1,  J ) I*BWM(I1,J)  *AWZP(I1,  Jl 
1)  -OWKP(Il.J) M AWM( II, Jll  +AWZ(I1,J1> ) -AWMII1, Jl) *BWZP(I1,J) ) 

20  CONTINUE  • • 

PEl=PEH-C3<;3U*1H22/2.*EU-ULl»£2/2.-C2*0L2»E3Ml.-XNl)*U0*C2*E‘»-E5 

PSYU=PSYU+0A4*C2*E6-DA3*£8 

pSAU=PSAU*E9/4. 

PSAG=C1*E10 
£ 1 - 9 • ' 

E2=0. 

E3=  0 . 

E4  = 0 . 

E5=  0 • 

£6=0. 

£8=0. 

00  21  J = 1 1 KK2 
JS= J**2 
JS  2 = J S*  * 2 

E1  = £10S2*  (CFM  ( 1 1 » J)  * * 2*-DFM  ( II  * J)  **2) 

E 2 =£2  +CF  MPP ( II ,J)**2+DFMPF ( 1 1 » J ) * *2 

£ 3=E J*JS*  ICF.H  II,  J)  *CFMPP(  11  , J)  4-0FM<  II,  J>  * UFMPP  ( 1 1 , J) ) 

E4=t*f JSMCFNPC I1,JI**2*DFKP(I1,J>**2) 

E5=£6*JB*CF1(I1 , J T 
tc  = E6*l,FfiPP(  II,  J) 

-_a  = £l»J*OFMP(Il*J) 

<?l  CONTINUE 

1 J*022/2.*E1*-D11/2.*£2-C2*0a2*E3*C2*E4/  ( ( 1 . -XN I ) * EXXP) 

- , rU-P  ;,VU*C2*UA  1*£5-DA2*E6 
• » YU»L1*2.  / U1.-XNI)*EXXP>  *E3 

T ' J ( ♦ P£  1 •£* 


’T*5 


STRYU=STkYU*PSY  U*E7 
SI  < AU=5  T R AU*-PS  A U *E  7 
3TRYG=STRYG*PSYG*£7 
STRAG  = STRAG«-PSA5*E7 
1,0  CONTINUE 

STRYU  = jA1*XNXX*-STRYU*OELTA/XL 

STRAU  = DAl*XNXX*STkAU*OEl_TA/XL 

ST  R Y G =2 . *XNX Y/ ( ' 1 . -X NI ) * EX XP ) -ST RYG* DELTA/ XL 

STRAG=2.  *XNXY/  ( ( 1 . - XN  I ) * E X XP ) - S T RAG*  [>E  L T A/  <2.*XL» 

P0T  = 3.1A159*RRMP0T*DELTA*XL*(U22*xNXX**2f2./((l.-XNI)*EXXP)* 
1XNXY’ *2-2.*STRAU*XNXX-2.*STRAG*XNXY) ) 

RETURN 

£NO 


SUBROUTINE  ADCG ( IEQ, Ml , CF , BF , AF , GF.NRHS, XNXX , XN XY ,LN, N J »NA W , 
lNJWtNF) 

C J 1.10N/30UND/LS  1 , LSN 

C Oil. 'ION/ r IDFR/QEL  TA,  AL1  ,GAl,AL2, 3T2,GA2 
GOMMON/FOURi  '<F O JR , K 1 , K2 , K 3 , K4 , K5 ,Kb , K7 , K8 , K9 
CO-1MON/FOUR2/K10,K11,K12.K13,K14,K15,KK2 
COilMON/CINTG/NEQPOT.MI  (500) 

DIMENSION  AF (Ml  * Ml ) . BF (Mi, Ml) ,CF(M1, Ml) ,GF (Ml. NRHS) 

C LSI  KIND  OF  BOUNDARY  CONDITION  OF  POTNT  1 

C LSN  KIND  OF  BOUNDARY  CONDITION  OF  POINT  NP 

C NAW  MAXIMUM  <*1  FOR  DIMENSION  OF  A 

C NBW  MAXIMUM  K FOR  DIMENSION  OF  b 
G NF  MAXIMUM  2*K  FOR  OIMEUSION  OF  C AND  D 

IF(IEQ.GT.l)  GO  TO  10 

CALL  RSTG(3F»CF, AF , GF , 1 , XNXX , XN XY , Ml , N J, NA W, N3 W ,NF , LN, NRHS) 

DO  2 11=1, K15  * • 

GF(I1+K15, NRHS) =GF(I1, NRHS) 

DO  2 J1=1,K15 

JF(Il*<i5,Jl)  = AJ_2*BF(Il,Jl)  + ALl*CF(Ii,Jl) 

DF (Il*K15,Jl*K15)=dT2*BF (II, J1)*AF(I1,J1> 

JF(  II  , J1*K15)=GA2*GF(I1,  Jl ) +GA1 *CF ( I 1 , Jl) 

2 CONTINUE 

CALL  dOUNDR(AF»CF»GF,1»XNXX,XNXY»LS1»M1*NJ ,N AW , NBW , NF , LN, NRHS ) 
DO  3 11=1, K15 
DO  3 J1=1,K15 
JF ( II , Jl ) = AL1*AF (II, Jl) 

AF (II +K15, j£)=dF( II , J1+K15) 

BF II1*J1*K15)=CF  (II, Jl) 

AF (II , J1)=GA1*AF(I1, Jl) 

3 CONTINUE 
RETURN 

10  IFdEQ.GT.2)  GO  TO  20  - 

CALL  RSTGIAF ,CF , BF , GF , 2 , XNXX , XNXY , Ml , N J , NA W , NBW , NF, LN, NRHS) 

DO  .4  11  = 1,  K15 
DO  4 J1=1,K15 

3F(I1  ,J1)  = BF  (II  , Jl)  f3T2*AF  (II,  Jl) 

GF(I1,J1*-K15)  = AL2*AF(I1,  Jl ) ♦■  AL1  *CF  ( 1 1 , Jl) 

AF  ( II  ,J1)=GA2*AF  ( 1 1 , Jl  > i-G  A1  *CF  ( II  ,J1) 

CF(I1,J1)=0. 

4 CONTINUE 
RETURN 

20  IF (ItQ.GE.NEGPOT-1)  GO  TO  30 
JP=IEQ 

CALL  RSTG(AF,CF  ,.3F,GF,  JP,  XNXX,  XNXY,  M1,NJ,  NAW,  NOW,  NF,LN,  NRHS) 

DO  5 11=1, K15 
DO  5 Jl  = l » K15 

3F(  II  , Jl  ) = DF  (II  , Jl)  +BT2*AF  (II,  Jl) 

TEMP=  GA  2 * A F ( I1,J1)«-GA1*CF(I1,J1) 

CF (II ,J1)=AL2*AF (II, Jl) f AL1*CF(I1,J1) 

AF  (II  , J 1 ) = T £ MP 

5 CONTINUE 
RETURN 

30  IF(IEO.EQ.NEQPOT)  GO  TO  40 


JP=IcQ 

. CALL  RSTGIAF.CF ,3F,GF, JP.XNXX,  XNXY,M1,NJ,NAH,N3W,NF  ,LN ,NRHS) 

JO  -6  11  = 1,  K15 
00  6 Jl  = 1 , K1 5 

.BF(Il,Jl)  = JF(Il,Jl)*ai2*AF(Il,Jll 
TEKP=GA2*AF( II, J1)*GA1*CF(I1,J1) 

CF(I1  ,J1)=AL2*AF  (11,  J1I»AL1*CF<U*J1I 
AF ( I 1 ,J1)=TEMP 
AF ( Ii, Jl+K15)=0 . 

6 CONTINUE 

RETURN  •- 

40  JP=IECl 

CALL  RSTG(AF,CF,3F,GF, JP , XNX X, XNX Y , Ml , NJ, N AW , N3 W, NF , LN , NRHS) 

00  7 11=1, K15 

GF (Il*K15,NRHS)=GF(Il»NRHS) 

JO  7 Jl  = 1 » K15 

3FII1 , Jl)=dF (II, J 1 » +BT2*AF (II, Jl> 

3F (II, J1*K15)=GA 2* AF(Il,Jl)*GAl*uF(Il,Jl) 

3FII1+K15,J1J=AL2*AF (II, Jl ) ♦ ALi *CF ( I 1 , Jl) 

7 CONTINUE 

CALL  BCUNDR(CF, AF , GF , J P, XNXX , XN XY , LS N, Ml , N J, NAW , N8W , NF , LN , NRHS) 
JO  d 11=1, K15 
TEMP=GF C II , NRHS) 

GF ( II , NRHS) = GF ( I1*K15»NRHS) 

GF (I1*K15»NRHS) = TEMP 
00  8 Jl= 1 » K15 

JF(IH-K15,J1*K15)=GA1*CF(I1,J1) 

CFCIH-K15, J1)=AL1*CF (II, Jl) 

CF  (II  ,J1)  = BF(I1+K15, Jl)  . . 

BF( II +K 15, J1)=AF (II, Jl) 

8 CONTINUE 
RETURN 
ENO 


SUBROUTINE  INVERT  (NA,  A ,C  ,M,Nf11 , NM 2 , Dt T , IX P , IDE T I 
UIHENSION  A INM1,NM1)  ,C  (NM2  > , M(NM2) 

DET=1. 

IXP  = G 
N.\'=NA 

IF (NN.NE.l)  GO  TO  303 
U£  T=  A ( 1 1 1) 

A(l,l)=l«/A(l,l) 

GO  TO  304 
3C3  00  9G  1=1, NN 
9 0 H C I > = - 1 

00  me  11=1, NN 
0=0 . D 0 

DO  112  K=1 , NN 

IFCM(K))  100,100,112  . 

100  00  110  L=1 , NN 

IF(MIL) ) 103,103,110 
103  IF (AbS(O) -AOS (A  K,L) ) > 105, 105,110 
10  5 LO=L  . . 

KO  = K 
0=  A ( < » L ) 

3 1 & A=  J 
110  CONTINUE 
112  CONTINUE 

IF(O.EQ.O.DO)  GO  TO  170 
GO  TO  138 

170  WRITE  16,502)  _ - 

STOP 

50  2 FORMA  T (/ , 5X, "DE TERMI N ANT  = 0 "/) 

188  NEMP=-M(LO) 

IKLO)  =M  (KOI 
M<KD).=  NEMP 
00  114  1=1, NN 
C(I)=A(I,LO) 

A(I,LD)=A(I,KO) 

114  A(I,KQ)=O.OQ 
A ( KD , KO ) =1 • OO 

• 00  115  J=1 , NN 

115  A IKU, J) =A  (KO,J) /O 
00  135  1=1, NN 
IF(I.EG.KO)  GO  TO  135 
DO  134  J=1 , NN 
TEHP=C(I>*A(KD, Jl 

134  A ( I»j)=A(I,J)-TEMP 

135  CONTINUE 

I F C IOET.NE.l)  GO  TO  140 

DET=UET*3IGA 

IF<KO.NE.LO)D£T=-OET 

629  IF  (AJS(DET)  .LT.  1.EHC)  GO  TO  630 
DET  = D!ZT/1.£«-10 

IXP=IXP*10 
GO  TO  629 

630  IF (ABS(DET)  .GT. 1.E-1&)  GO  TO  140 


CET  = DET*l.t»10 
. IXP=IXP-10 
140  CONTINUE 

JO  200  I-l.NN 
. L=0 

15  C L = k.*  1 

IF(fKL)-I)  150,  160,150 
160  F(l)=M(I) 

,1(11-1 

JO  200  J=1,NN 
TEMP=  A (L  , J 1 
A (L, J)=A(I, J) 

200  A(I,J)  = rEI1P 
30  4 RETURN 
£NO 
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SU3R0UHNE  YMY  (Nl,  A , B , C,  N2 , 1 1 , 2,  L3 , T) 
DIMENSION  A(L1,L2I,B(L1,L1),C{L1,L2),T(L3» 
IF (N2  .cQ.ll  GO  TO  100 
00  11  1=1. N1 
00  10  J = 1 » H 2 . 

T E II  ? = 0 • 


00  20  K=1,N1 

20  TEN P=TEMP*StItK)  *C ( K*  J ) 

10  I (J)=TEMP 

00  3u  J=1*N2 
30  A ( I . J ) =.T  ( J ) 

11  CONTINUE 
RETURN 

100  00  111  1=1, N1 

T£M?= 0 • 

00  120  K=1,N1 

120  TEHP=T£:-lP«-3(I,K)*C(Kf  1) 
111  TII)=TEMP 

00  130  1=1, N1 
130  A(I,1I  = T II) 

RETURN 
END 


r 
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SUURU  UT I N£  YSYMY  1N2 .N 1 , A ,B,C . D , N3 , Ll , L2 . L 3 , L4 ♦ II 
DIMENSION  A(L1, L 3) ,9(L1,L3) ,C(L1,L2> ,D (L2.L3) , T (L4> 
IF (N3.EQ.1)  GO  TO  1G0 
DO  11  1=1, N1 
DO  10  J = 1 , N 3 
TEMP=0. 

DO  20  K=  1 » N2 

20  TErtP=TEMP*CtI,K) *DIK,J> 

10  T(J)=or4I,J)  - TEMP 
DO  30  J = 1 » N3 

30  A(I,J)=T(J> 

11  CONTINUE 
RETURN 

100  OO  111  1=1, Ni 


TEMP  = 0 • 

OO  120  K=1,N2 

120  TE.3P=TEMP+C(I,K)*D(K,1) 
111  T(I)=0(I,1>-TEMP 
OO  130  1=1, Nl 
130  A ( I , 1 ) =T  II) 

RETURN 

END 


i 

i 
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SUBROUTINE  XREAO(Nd, A , LI , L2, Ml , M2 , INO, M3 , V V) 

COMMON/COISK/I2 1(501) ,122(501) 

RECORD^IND  OFMOIREcfVACCESS  OATA  SET  NO  IS  READ  AND  ALLOCATED 
JY  ROWS  INTO  L1*L2  PORTION  OF  MATRIX  A 
L3=L1*L2 

CALL  READMS(ND,Vi/»L3»IND)r .....  . 

<L  - 0 


DO  10  NROW=l ,L1 
30  10  NCOL=t,L2 
' X L = KL+1 

A(NROW,NCOL)=VV (KL) 
10  CONTINUE 
RETURN 
END 


fc3 


SUOROU1 1NE  XWRI TE( NO. A, LI. L2, Ml  , M2, 1 NO, M3, Vtf> 

• C01M0N/CDISK/I2K501)  , I22CBQ1) 

DIMENSION  A(M1 , M2) ,VV (MS)  , _ , 

L1*L2  PORTION  OF  MATRIX  A IS  WRITTEN.  BY  ROWS  ON  OIRECT  ACCESS 
DATA  Scl  NO  IN  RECORD  IND 
<L-  0 

00  1J  NR0W=1,L1 
00  10  NCCL=1,L2 
<L  = KL  *1 

VV(KL)-A(NROW'»NCOLI 
10  CONTINUE 

CALL  WRI  TI1S  (NO,  VV  ,Kl,  INO.-l) 

RETURN 

END 
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SJJkOUJ  IN£  POT£kS( IOt  T ,NkHS,MAXN, AP,bP,CP, GP,PR.XP,C,N  T,T1, 

1 Vl,lA<2,IXPM,  D£ TM.XNXX, XNXY,LN,NJ,NAW  ,NBW,NF> 

C MAX2=MAXN*MAXN 

COM t'. ON/C  INF G/NE  QPO  T , Ml  15  00) 

COMMON/ GUISK/I 21 (501) ,122(501) 

DIMENSION  AP(MAXN.MAXN)  , BP  (MAXN  , MAXN)  , CP  C M AXN , M AX  N) 

DIMENSION  PR  (MAXN,  MAXN) , GP (MAXN , NRHS > , XP(MAXN, NRHS ) 

DIMENSION  T1 (MAXN) ,C (MAXN) , Ml (MAXN) , VI (NAX2) 

C EQUIVALENCE  (AP ( 1 , 1 ) , V 1 ( 1) ) 

I XPM=  0 
QETK=1. 

00  100  I-l,N£QPOr 

CALL  AbCG(I,MAXNtCP,BP,AP,GP,NRHS,XNXX,XNXY,LN,NJtNANf N8W.NFI 
N=  M I ( 1 1 

IF(I.cQ.l)  GO  TO  88» 

NMINl-=MI  (1-1) 

888  IFd.cQ.NSGPOT)  GO  TO  999 
NPLUS1=M I (1+1) 

999  CONTINUE 

IF(I.cQ.l)  GO  TO  12 

CALL  YSYMY { NMI  N 1 , N , BP , BP , CP , PR , N, MAX N, MAXN ,MAXN , M A XN , T 1 ) 

12  CALL  INVERKN,  BP,  C,MT,  MAXN,  MAXN,  QET , I XP , I OtTI 
IF(IOET.NE.l)  GO  TO  640 
Oc  T M = Uc  T * D E T M 
IXPM= IXP+IXPM 

IF(ABSIDETM)  .LT  .l.E+10)  GO  TO  630 
l)£  TM=  Dt  T M/ 1 • E *-1  0 
IXPM=IXPMHO 
«0  TO  640 

630  IF ( A3S (DETM) .GT , l.E-10)  GO  TO  640 
OcTM=DETM*l.E*l  0 
I XPM= I XPM-1 0 
640  CONTINUE 

IF ( I . EQ . NEQPOT)  GO  TO  102 

CALL  YMY(N, PR,  BP, AP,NPLUS1, MAXN, MAXN, MAXN, Tl) 

CALL  XWRIT£(21,  PR,N,NPLUS1 ,MAXN  , MAXN , I , MA X 2, VI) 

1C2  IF ( I . £Q. 1)  GO  TO  32 

CALL  YSYMY (NMIN1»N,XP, GP,CP,XP, NRHS , MAX N, M AXN, N RH S , MAX N , T 1 ) 
CALL  YMY l N , XP,  BP , XP , NRHS , M AX N , N RHS , M AXN , T 1 ) 

GO  TO  42 

32  CALL  YMY (M,XP,3P,GP, NRHS, MAXN, NRHS, MAXN, Tl) 

42  CALL  XNRITE(22, X?,N, NRHS, MAXN, NRHS, I, MAXN,  Tl) 

100  CONTINUE 

ME  QPO  T = NEQPOT- 1 
00  200  K=l,MEQPOr 
NK  = N£QPO  T-K 
NMIN1=MI (NK) 

N=MI (NK+1) 

CALL  XRE AO (21,  PR, NMI Ni ,N, MAXN, MAXN, NK, MAX 2, VI) 

CALL  XRE AD (22,  GP,NM INI , NRHS, MAXN, NRHS, NK, MAXN, T 1) 

CALL  YSYMY (N,NM INI, XP.GP, PR, XP , NRHS, MAXN, MAXN,  NRHS, MAX N,T1) 
CALL  XWR1TE(22, XP,NM INI, NRHS, MAXN, NRHS, NK, MAXN, Tl) 

200  CONTINUE 


SUBROUTINE  TkANSF(Tl,MAXN, IDER, IPRR) 
•COMMON/FOURl/KFOUR,Klf  K2,K3,K4, K5  ,K6 , K 7 , K 8 ,K9 
C0M'10N/F0UR2/K1U,K11,K12  .K13.K14 ,K15 ,KK2 
COHMON/OINTG/NEUPOT ,HI  (500) 

SOMN ON/GUI SK/I21 (501), 122(501) 
C0MMON/FIOFR/OELTA,Al1,GA1,AL2, BT2,GA2 
COMrlON/PkCSl/AWI-HlOO,  3)  , AW MP  ( 1 0 0 , 3) , AW  MPP  (10  0,3) 
COMMUN/PKES2/8WM ( 10  0 , 2 ) , SW MP  ( 1 0 0 » 2)  , 3WMPP (100,2) 
C0I1M0N/PRES3/CFM  (10  0,4)  ,CFMP  (10  0, 4) , CFMPP ( 10 0 , 4 ) 
COilMON/PxE  S4/13F  M(1G0*4),  UFMP  ( 1 0 0 , 4)  , DF MPP  (100,4) 
COMMGN/PRES5/AW  Z(100,3)» AW ZP ( 1 0 0, 3) , AWZPP (100,3) 
COMMON/PRE S6/3WZ (100, 2)  , 3WZP  (1  u 0,  2) , 3WZPP(  10  0,2) 
C0MM0N/RES9N1/A WMB(2, 3), AWMPP3(2, 3) , uWMB(2 ,2) , BWMPP£3(2 ,2) 
COMMON/R£S3N2/CFMB(2,4  ) , CF MPPB ( 2, 4 > » Or  MB ( 2 , 4 ) ,OFMPPB(2 ,4). 

OI. HEN  SION  f 1 (HA  XN)  

IF ( IPRR. EQ. 3)  GO  10  278 
00  10  Il=l,NEQPOr. 

NL=MI(I1) 

. CALL  RE A QMS ( 22 , Tl,NL»Il>  . _ . 

I F ( 1 1 . N £ . 1 ) GO  TO  175 
30  11  Jl=l ,K1 
AWM3(1, Jl)=Tl( Jl) 

AWM(I1,  Jl)  = ri(JlfK15) 

AWMPPB ( 1 , J 1) =T1 (J1  + K7) 

AW  MPP  (I1,J1)-=T1  ( J1+K7*K1 5) 

11  CONTINUE 

00  12  J1=1,KF0UR 
3WH31 1, Jl)=Tl( Jl+Kl) 

BWM(I1,  J1)=T1(  JH-K1H<15) 

8WMPP9 ( 1 , Jl) =T1 ( J1+K9) 

3W MPP (II ,J1) =Tl ( J1+K9+K15) 

12  CONTINUE  ' ~ 

00  13  J 1=1 , KK2 

CFH3(1,  J1)=T1(J1+K3) 

OFHB(  1,  J1)=T1(J1K5) 

CFM(I1,J1)=T1(J1+K3*<15) 

OFMIII,  J1)=T1(  J1  + K5*K15) 

CF.MPPB  ( 1 , J 1)  =T1  (J1*K11) 

0FHPPB(1,J1)=T1  (JH-K13) 

CFMPP  (I1»J1)=T1  (J1*K11*K15)  ^ 

0FMPP(I1,J1)=T1(J1  + K13«-K15) 

13  CONTINUE 
GO  TO  10 

175  DO  14  J 1 = 1 ,K1 

AWM(I1.  J1)=T1(J1) 

AWMPP ( II  ,J1)=T1  (J1+K7) 

14  CONTINUE 

00  15  J 1 -1 , KFOU R 
BWMdl,  J1)=T1(J1«-K1) 

3WMPP (II , Jl) =T1'(  J1  + K9) 

15  CONTINUE 

00  16  J1=1,KK2 
CFM(I1,  J1)  = T1(  J1*K3» 
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JFM(I1, J1)  = T1(J1*K5) 

C F M ?P  (IltJl)-Tl  (J1*-K11) 

OF  MPP (II , J 1 ) -T 1 ( J1  + K13) 

16  CONTINUE 

IF ( II . NE .NEOPOT ) GO  TO  10 
DO  17  J1=1,K1 
AWM3(2,  Jl)  = T1(  J1*K15) 

AWMPPB(2 * J 1 ) = T 1 ( J1+K7*K15> 

17  CONTINUE 

QO  18  Jl  = l, KFOUR 
dW?1G(2,  Jl)  =T1(  JH-K1K15) 

GWMPP3(2  ,J1)=T1  ( J1«-K9*K15) 

18  CONTINUE 

00  19  J1=1,KK2 
CFNB(2,J1)=T1( J1*K3*K15) 

OFMC(  2,J.1)=T1  (J  i*K5«-K15) 

CFNPPB<2 ,J1)=T1 (J1+K11+K15) 

DFMPP  3(2  ,J1)  =T1  ( J1K13+K15) 

19  CONTINUE 
10  CONTINUE 

IF (IDcR.NE.l)  GO  TO  275 
NEUP=NEJP0T-1 
00  20  11=2, NEQP 
' oo  21  J i = l , :<  i 

Afc.IPl  II,  J1)  = AL1*AWM(I1-1,J1)  +GA1*AWM  (Il  + l , JD 
IFCJ1.EQ.K1)  GO  TO  21 

8WMP ( II, J1)=AL1*BWM( 1 1-1, Jl) +G A 1*  3WM ( 1 1 >1 ♦ Jl) 

21  CONTINUE 

00  22  Ji= 1 , <K2 

CF1P  ( Il»Jl)  = ALl*CFM(Il-ltJl)  +GAl*CFN(IlHtJl) 
DFMPdl,  Ji)=ALi*DFMdl-l.Ji)  *,GA1*0FM(I1  + 1»J1I 

22  CONTINUE 

20  CONTINUE 

00  23  Jl  = l ,'K1 

AWMPl 1, Jl) = AL1* AWMUd, Jl ) ♦ G A 1* A WM  ( 2 , Jl) 
AWnPTNEQPOT.Jl)  = AL1*  A WN  ( NEQP  , J 1 ) «-GAl  * A WM8  ( 2 , Jl ) 
IF(Jl.EQ.Kl)  GO  TO  23 

BWiTPd,  Jl)-ALl*0rlM5(l»  Jl ) fGA 1 * BWM  (2,  Jl) 

BWMP  ( NEQPO  T , Jl)  = AL 1* 3W M ( NEuP  , J 1 ) «■  GA1* B WHB  ( 2 , Jl ) 

23  CONTINUE 

DO  24  J1=1,KK2 

CFMP(1 , Jl) =AL1*CFMB ( 1, Jl ;+GAl*CFM (2, Jl) 

DFMP( 1, Jl) =AL1*  DFMBd,  J1)«-GA1*LFH(2,  Jl) 

CFMP ( NEQPO T , Jl ) = Al1*CFM ( NEUP , J 1 ) >GA1*CFMB ( 2, Jl ) 
OFMP  ( NEOPOT  , Jl ) = AL1*DFM(NEQP, Jl)*GAl*DFHd (2, Jl) 

24  CONTINUE 

275  IFdPRR.NE,  1 ) RETURN 
278  CONTINUE 
J1  = 0 

WRITE (6,400) Jl 
WRITE (6,500) 

X X = 0 . 

WRITE (6, 60  0) AW  MB (1, 1) , AWMPP3 (1,1) 


00  48  11  = 1,  NEUPOT 

WR{  Tfc  (6,509)11,  XX,  AWM(I1,1*  , AWMP(I1,1)  ,AWMPP(I1,1I 
XX=XX  +GELTA 
<♦8  CONTINUE 

. WRIT-  (6,600) AWMB(2,1)  , AWMPPB(2,1) 

DO  49  J 1 = 1 , KK2 
IF (Jl.GT  .KFOUR)  GO  TO  68 
WRITE  (6, 400  Ji 
WRIT! (b, 500) 

WRITE (6.7G0) AWM311, Jl*l)  .AWMPPB  (1  , Jl  +1 ) , i3W  M3  (1 . Jl ) , 3WMPP3  C 1 , Jll 
00  51  I1=1,N£QR0T 

WRITE  (6, 609)  11,  AWM(I1,  Jl  + 1)  , AWMP(  II,  Jld)  , AWMPP (II , Jl* 1* . 
laWMIIl, Jl) » 3WMP ( II , Jl) ,3WNPP(I1 , Jl) 

51  CONTINUE 

WRITE  (6,700)  AWMBI2,  Jl*l>  , A WMPP3  (2  , Jl *-1 ) , 3W  MO( 2 , Jl ) , 8WHPPB (2 , Jl > 

6 8 WRITE  (6,4 0 0) Jl 
WRITE (6,501) 

WRITE  (6, 7 00)CFMB(1, Jl)  , CF MPP B ( 1 , J 1J , UFMO( 1 ,Jl) , DFHPP3(1,J1) 

DO  52  11  = 1 , NEQP  OT 

WRITE (6,60  9) 11 , CFM ( 1 1 , Jl ) , CFMP ( II , Jl ) , CFMPP( 11 , Jl) , OFM ( II , Jl ) , 
1DFMP(  II,  Jl) , DFHPPdl,  Jl) 

52  CONTINUE 

WRITE(6,70  0)CFMB(2,J1)  ,CFMPP3(2  , J 1 ) , OF M3 ( 2 , Jl ) , OF HPPB ( 2 , Jl ) 

49  CONTINUE 

40  0 FORMAT!//, 2X, "RESULTS  FOR  KFOUR=",  I8/2X,  "***•**•*•**+**•*•**•*+• 

^****•**44**4444*  444444444  •«} 

50  0 FORMAT  (//,2X,"P 0INT",4X,  "LENGTH”,  9X , "WCOS”  , 11X  , "WPCOSM , 10X  , 

1"W PPG  OS",  10X,”WS IN”, 1 IX, "WPS IN”, 9X,”WPPSIN"/2X, ”==  = ==  ======  = ===== 

2===== ========================== "7 

501  FORMAT ( //, 2X ,"P  0 IN  T" , 4 X , "LENGTH", 9X, "FCOS" , 1 IX , "FPCOS" , 10 X, 

1”F  PPG  OS",  1 OX,  "F  SIN",  1 IX,  •JFPSIN”,  9X,"FFPSIN*V2X  ,"==  ===  = = ==  = =====  = = 

2====================;===========”) 

600  FORMAT (//20H  FICIIVE  POINT  E 1 5 . 6 , 15X ,E 15 . 6//) 

5C9  FORMAT  (18, El  2.  4, 3E  15,6) 

700  FORMAT (//20H  FICTIVE  POINT  E15 .6 , 15X ,2£1 5 . 6, 15X, E15. b//) 

609  FORMAT  (I6,12X,6E15.6) 

RETURN 

END 
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SUBROUTINE  IMPERF 

COMMON/FOURl/KFOUR,Kl, K2,K3,K4, K5,K6,K7,K8,K9 
C0MM0N/F0UR2/K1 0 ,K1 1 ♦ K12 , ML  3 , K1 4 , K15 , KK2 
COMMON/C INT3/NE  QPOT » Ml (5  00) 

G0MM0N/FACT3/XL.XH 

COMMON/F IOFR/ DELTA, AL1»GA1,AL2»BT2»GA2 
CO M MON/P RE S5/AWZ(10O»3)»AWZP (100,3), AW  ZPP ( 10  0*31 
C0MM0N/PRES6/BWZ (100,  2), 3WZP (10  0, 21 , 3WZPP (10  0,2) 
DIMENSION  AM (10) , BM ( 10 ) 

DV  = 1. 

PI  = «*.*ATAN(DV> 

READ(5,20P)ML,DL 

RE  AD  (5, 100 1 (All  (I  I,  1=1,  ML) 

READ( 5,100) (DM (I), 1=1, ML) 

OO  12  1=1, ML 
AM(I)=AM(I)*OL*XH 
BM(I)=BM(I)*XH*DL 
12  CONTINUE 
200  FORMAT (16, E12. 4) 

100  FORMAT (6E12. 4) 

XX=0. 

UO  1C  11  = 1 » NEQP  OT 
AWZ  (11,21  = 0. 

AWZP( Il,2)=0. 

AWZPPTIl ,2)=0. 

8WZ(11,1)=0. 

3WZP(I1,1)=0. 

awzpp(ii,i)=o. 

AMZ(I1,1)=0. 

AWZ P ( ll,l)=0. 

AW  ZPP  (I1,1)=0* 

BWZ(I1,2)=0. 

OWZP(I1,2)=0. 

3WZPP(I1,2)=0. 

DO  11  IM=1 »ML 
A1=PI*IM/XL 
A 12  = A 1**2 
A1X=SIN(A1*XX) 

A2X=C0S( A1*XX> 

AWZ(I1,2)=AWZ(I1,2) -AM (IN) *A1X 
AWZ?(I1,2)=AWZP(I1,2)-AM(IM)*A1*A2X 
AWZPP  (I1,2)=AWZPP(I1,2)*AM(IM)*A12*A1X 
3WZ(I1,1)  = 3WZ(  I1,1)-3M(IM)*A1X 
‘ BWZP( II, 1) =3WZP (I1,1)-BM (IM) *A1*A2X 
dWZPPCIl,l)=BWZP?(Il,l)*bM(IM)*Al2*AlX 
11  CONTINUE 

IFIK1.LE.2)  GO  TO  50 
DO  13  Jl  = 3 »K1* 

• A W Z ( 1 1 , J 1)  = 0 • 

AWZPI II » J1)=0, 

AWZPP(I1,J1)=0. 

3WZ(I1,J1)=0. 

BWZP( II, J1)=0. 


map*** 


UWZP*MI1,J1)=0 
CONTINUE 
CONTINUE 
XX=XX«-DELTA 
CONI INUE 
*ETU<*N 
ENO 


APPENDIX  B 

Summary  of  Work  Performed  Under  AFOSR- 74-2655 


The  overall  effort  under  grant  AFOSR- 74-2655  which  deals  with  the 
minimum  weight  design  of  fuselage- type  stiffened  cylindrical  shells,  can 
best  be  described  by  giving  a precise  but  general  statement  of  the  pro- 
blem and  by  discussing  individually  all  of  the  accomplishments. 

The  precise  statement  is  as  follows:  given  an  internally  stiffened, 
imperfect,  circular,  cylindrical  shell  of  specified  material,  radius,  and 
length,  find  the  size,  shape  and  spacings  of  the  stiffeners  and  the  thick- 
ness of  the  skin,  such  that  the  resulting  configuration  can  safely  carry 
a given  set  of  destabilizing  loads  (applied  individually  or  in  combination) 
with  minimum  weight.  The  solution  of  the  problem  required  application  of 
modern  optimization  techniques,  development  of  efficient  solution  methodo- 
logies for  finding  critical  conditions  for  stiffened  shells  under  certain 
loads  and  in  the  presence  of  initial  imperfections,  and  incorporation  of 
all  of  the  above  into  a single  computer  code. 

The  load  cases  considered  are  (a)  uniform  axial  compression,  (b)  uni- 
form or  nonuniform  lateral  pressure  and  (c)  uniform  torsion.  These  loads 
are  representative  of  the  types  that  the  configuration  is  expected  to 
encounter  in  service  when  used  either  in  an  aerospace  vehicle  fuselage  or 
in  a submarine  hull. 

The  structural  theories  employed  in  the  mathematical  model  are  based 
on  the  following  assumptions. 

(1)  the  material  behavior  is  linearly  elastic. 

(2)  the  kinematic  relations  correspond  to  those  of  moderate  rotations 
[(rotation)^  « l] 

(3)  the  loads  are  applied  quasi-statically 

(4)  the  stiffener  spacings  are  small  and  the  connections  are  monoli- 
thic so  that  the  "smeared"  technique  be  applicable 

(5)  the  imperfection  shapes  are  smooth  functions  of  position  and  the 

imperfection  amplitude  is  reasonable  from  a manufacturing  point 
of  view  (several  skin  thicknesses). 


The  overall  effort  can  be  divided  in  two  parts  (a)  those  problems  which 
are  free  of  initial  imperfections  and  for  which  linear  buckling  analyses 
have  been  employed  in  the  optimization  technique  and  (b)  those  problems 
which  include  geometric  imperfections. 

The  accomplishments  for  each  class  of  problems  are  listed  separately. 

OPTIMIZATION  OF  PERFECT  GEOMETRY  CYLINDERS 
EMPLOYING  LINEAR  BUCKLING  ANALYSES 

The  work  associated  with  this  phase  of  the  research  program  has  been 
reported  in  detail,  through  three  AFOSR  Technical  Reports  (Refs.  B.1-B.3) 
and  numerous  publications  in  refereed  journals  (Refs.  B.4-B.11). 

In  this  phase,  the  design  objective  of  the  aforementioned  optimization 
problem  is  minimum  weight . The  general  instability  load(s)  is  considered 
to  be  the  one  or  one  of  the  active  failure  modes  (in  some  cases,  applicable 
to  pressure  loaded  submarine  hulls,  skin  yielding  is  considered  to  be  the 
active  failure  mode).  Thus,  in  the  optimization  formulation  this  load  is 
taken  as  an  equality  constraint.  All  other  failure  modes  are  taken  as 
inequality  constraints  (behavioral).  These  include  panel  instability,  lo- 
cal instabilities  of  the  skin  and  stiffeners  and  yielding  of  the  skin  and 
stiffeners.  Other  inequality  constraints  are  of  the  geometric  type  and 
they  represent  realistic  dimensions  for  some  of  the  design  variables  such 
as  minimum  gages  for  skin  and  stiffener  thicknesses  as  well  as  limitations 
on  the  stiffener  spacings. 

Moreover,  because  of  the  possibility  of  detrimental  effects  due  to 
failure  mode  interaction,  the  condition  of  separation  of  these  modes  is 
imposed  on  the  solution. 
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The  design  variables  of  the  problem  are  the  skin  thickness,  the 


stiffener  spacings , the  stiffener  geometries  (shapes  considered  are: 
rectangular,  R,  angle.  A,  T-shapes,  T,  channel  sections,  C,  I-sections,  I, 
hat-sections,  H,  and  others),  web  and  flange  widths  and  thicknesses.  The 
number  of  these  variables  is  eleven^but  due  to  the  small  effect  of  the 
stiffener  web  and  flange  thicknesses  on  the  minimum  weight, it  is  assumed 
that  these  thicknesses  are  equal  and  thus  the  number  is  reduced  to  ten. 

The  load  cases  considered  and  reported  in  Refs.  B.l-B.ll  include  (i)  uni- 
form axial  compression,  ( ii ) uniform  pressure,  (iii)  torsion,  and  (iv)  all 
possible  combinations  of  above. 

Among  the  most  important  conclusions  of  this  phase  one  may  list:  (a) 
the  developed  methodology  includes  the  following  desirable  features: 

(i)  The  designer  can  easily  assess  the  need  or  lack  of  need  for 
stiffening  in  both  directions. 

(ii)  through  no  penalty  or  minimum  penalty  in  weight  the  designer  may 
avoid  failure  mode  interaction. 

(iii)  the  designer  may  carry  out  important  trade-off  studies  to  arrive 
at  a practical  minimum  weight  configuration. 

(b)  For  axial  compression  and  pressure  loaded  systems,  the  minimum  weight 
design  is  not  unique.  This  means  that,  for  a given  set  of  the  specified 
parameters,  the  design  variables  can  be  adjusted  so  as  to  give  several 
acceptable  designs  corresponding  to  the  same  minimum  weight. 

(c)  The  optimum  distribution  of  material  is  load  case  dependent  and  length 
to  radius  ratio  dependent.  For  moderate  length  cylinders  (l/R~  3),  the 
following  observations  are  made; 

(i)  The  optimum  distribution  of  material  corresponds  to  approximately 
60%  in  the  skin,  30%  in  the  stringers  and  10%  in  the  rings  for 
axially  loaded  shells. 

(ii)  For  pressure  loaded  shell,  the  optimum  distribution  of  material 
corresponds  to  60%  in  the  skin,  30%  in  the  rings  and  10%  in  the 
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s trlngers . 

(iii)  For  combined  axial  compression  and  pressure  the  combined  stiffe- 
ner material  is  approximately  407„  of  the  total,  but  distributed 
into  stringer  and  ring  material  in  accordance  with  the  relative 
amounts  of  pressure  and  axial  compression  (with  respect  to  li- 
near theory  critical  loads) . 

(iv)  For  pressure  loaded  shells  the  amount  of  ring  material  increases 
with  the  L/R  ratio  (slightly). 

(d)  The  optimum  stiffener  shape  is  also  load  case  dependent. 

(i)  For  axially  loaded  stiffened  shells  the  stringers  must  be  T- 
shaped,  while  the  rings  must  be  rectangular  (see  Fig.  B.1-B.3) 

(ii)  For  pressure  loaded  stiffened  shells  the  rings  must  be  T- 
shaped,  while  the  stringers  must  be  rectangular. 

(iii)  For  torsion  loaded  stiffened  cylinders  the  combination  of  hat 
stringers  with  either  hat  rings  or  rectangular  rings  proved  to 
be  the  most  efficient  one. 

(e)  The  curve  for  determining  the  optimum  skin  thickness  (part  of  the 
minimum  weight  design  methodology)  is  relatively  flat  for  the  axial  com- 
pression and  uniform  pressure  cases.  Therefore,  very  precise  determination 
of  the  optimum  thickness  is  not  necessary  for  minimum  weight  design.  In 
the  case  of  torsion,  the  minimum  weight  thickness  is  always  equal,  to  the 
minimum  gage  imposed. 

OPTIMIZATION  OF  IMPERFECT  STIFFENED  CYLINDERS 

The  work  associated  with  this  phase  of  the  program  has  been  reported 
through  Refs.  B.12-B.17,  as  well  as  in  the  present  AFOSR  Technical  Report. 

In  this  phase  a methodology  was  first  developed  and  demonstrated 
through  numerous  examples  for  analyzing  imperfect  stiffened  shell  configu- 
rations under  various  load  conditions  (see  Refs.  B.12-B.15  and  present 
report) . 

Then,  this  methodology  was  incorporated  with  the  linear  buckling  theory 
optimization  procedure  in  order  to  optimize  an  imperfect  stiffened  cylinder 
(see  Refs.  B.12,  B.16,  B.17  and  present  report). 
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Thus,  a solution  methodology  has  been  developed  in  order  to  optimize 
an  imperfect  stiffened  circular  cylindrical  shell  under  individual  and 
combined  application  of  destabilizing  loads. 

Finally,  a small  effort  was  exerted,  under  the  present  grant,  to 
investigate  the  imperfection  (load  eccentricity)  sensitivity  of  a dif- 
ferent structural  configuration  (a  simple  two-bar  frame--see  Ref.  B.18). 
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FigB.2Case  1-  Effect  of  Ring  Shapes  on  Cylinder  Weight 
using  Most  Efficient  Stringer  (TS  or  IAS.  Cx- 1X07 ) 
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